{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data downloading into ASDF and computing cross correlation functions \n",
    "\n",
    "In this notebook, we show the basic steps of NoisePy for 1) downloading seismic data using [ObsPy](https://github.com/obspy/obspy/wiki) functions (`get_station` and `get_waveform`) and saving them into ASDF file, and 2) computing cross correlation functions. As you will find through the practise, apart from the data downloading and data loading processes (sections 0-3), the rest of the notebook is essentionally the same to the one shown in the other notebook of `cross_correlation_from_sac`. If you would like to use the NoisePy to deal with SAC/mSEED files at your local disk or any other format that can be ready by ObsPy, we recommend you look at the other notebook first. \n",
    "\n",
    "The steps to compuate cross correlation functions as illustrated in this notebook are:  \n",
    "* Download Z-component data from 2 TA stations, pre-process the data (remove gaps, downsample etc) and save them into one ASDF file\n",
    "* Read the ASDF file with the first station to be source and the other as the receiver\n",
    "* Break the continous data into small segments with overlaps\n",
    "* Perform Fourier Transform to convert signals into frequency-domain\n",
    "* Calculate cross correlation functions between the small time segments and choose to stack (substack) the cross correlation function of each segment and return to time domain\n",
    "* Save the correlation function into an ASDF file\n",
    "\n",
    "More details on the descriptions of data processing, parameters for different cross correlation method and performance of NoisePy can be found in the online [documentations](https://noise-python.readthedocs.io/en/latest/) and the manuscript.\n",
    "\n",
    "`Jiang, C. and Denolle, M. NoisePy: a new high-performance python tool for seismic ambient noise seismology. In prep for _Seismological Research Letter_.`\n",
    "\n",
    "\n",
    "\n",
    "Chengxin Jiang & Marine Denolle\n",
    "\n",
    "Department of Earth and Planetary Science\n",
    "\n",
    "Harvard University\n",
    "\n",
    "November 2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building env for NoisePy\n",
    "\n",
    "Before running this notebook, make sure that you have created and activated the conda env made for NoisePy. If not, you can create one using command lines below (note that jupyter is installed with the command lines here in order to run this notebook). \n",
    "\n",
    "```python\n",
    "conda create -n noisepy -c conda-forge python=3.7.3 numpy=1.16.2 numba pandas pycwt jupyter mpi4py=3.0.1\n",
    "conda activate noisepy\n",
    "pip install obspy pyasdf \n",
    "git clone https://github.com/mdenolle/NoisePy.git\n",
    "```\n",
    "\n",
    "Then you need to activate this notebook with the newly built NoisePy env by invoking the jupyter with the following command line.\n",
    "\n",
    "```python\n",
    "jupyter notebook\n",
    "```\n",
    "\n",
    "Now we can begin to load the modules needed for downloading part first. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import time\n",
    "import obspy\n",
    "import pyasdf\n",
    "import os, glob\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import noise_module\n",
    "from obspy import UTCDateTime\n",
    "import matplotlib.pyplot as plt\n",
    "from obspy.clients.fdsn import Client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Settings for data request\n",
    "\n",
    "Here we download one day of continous data (Z component) for 2 TA stations located west of the Boston city, and we do not remove instrument response (set `rm_resp` to `inv` if you want to use the inventory for instrumental removal), downsample it to 1 Hz and filter it to 0.05-0.4 Hz in the pre-processing. Note that in NoisePy, the way of specifying network, station and channel are slightly different and here just show one example. NoisePy also accepts using a compliated station list for data downloading in S1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download parameters\n",
    "client    = Client('IRIS')                                      # client/data center. see https://docs.obspy.org/packages/obspy.clients.fdsn.html for a list\n",
    "samp_freq = 1                                                   # targeted sampling rate at X samples per seconds \n",
    "rm_resp   = 'no'                                                # select 'no' to not remove response and use 'inv','spectrum','RESP', or 'polozeros' to remove response\n",
    "respdir   = './'                                                # directory where resp files are located (required if rm_resp is neither 'no' nor 'inv')\n",
    "freqmin   = 0.05                                                # pre filtering frequency bandwidth\n",
    "freqmax   = 0.4                                                   # note this cannot exceed Nquist freq                         \n",
    "\n",
    "chan = ['BHZ','BHZ']                                            # channel for each station\n",
    "net  = ['TA','TA']                                              # network for each station \n",
    "sta  = ['K62A','K63A']                                          # station (using a station list is way either compared to specifying stations one by one)\n",
    "start_date = [\"2014_01_01_0_0_0\"]                               # start date of download\n",
    "end_date   = [\"2014_01_02_0_0_0\"]                               # end date of download\n",
    "inc_hours  = 24                                                 # length of data for each request (in hour)\n",
    "nsta       = len(sta)\n",
    "\n",
    "# save prepro parameters into a dic\n",
    "prepro_para = {'rm_resp':rm_resp,'respdir':respdir,'freqmin':freqmin,'freqmax':freqmax,'samp_freq':samp_freq,'start_date':\\\n",
    "    start_date,'end_date':end_date,'inc_hours':inc_hours}\n",
    "\n",
    "# convert time info to UTC format\n",
    "starttime = obspy.UTCDateTime(start_date[0])       \n",
    "endtime   = obspy.UTCDateTime(end_date[0])\n",
    "\n",
    "# another format of time info needed for get_station and get_waveform\n",
    "s1,s2 = noise_module.get_event_list(start_date[0],end_date[0],inc_hours)\n",
    "date_info = {'starttime':starttime,'endtime':endtime} \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Download data and save into ASDF file\n",
    "\n",
    "Note that this step might take a few seconds for data download. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "working on station K62A\n",
      "working on station K63A\n"
     ]
    }
   ],
   "source": [
    "# write into ASDF file: using start and end time as file name\n",
    "ff=os.path.join('./',s1+'T'+s2+'.h5')\n",
    "with pyasdf.ASDFDataSet(ff,mpi=False,compression=\"gzip-3\",mode='w') as ds:\n",
    "\n",
    "    # loop through each station\n",
    "    for ista in range(nsta):\n",
    "\n",
    "        # get inventory for each station\n",
    "        try:\n",
    "            sta_inv = client.get_stations(network=net[ista],station=sta[ista],\\\n",
    "                location='*',starttime=s1,endtime=s2,level=\"response\")\n",
    "        except Exception as e:\n",
    "            print(e);continue\n",
    "\n",
    "        # add the inventory into ASDF        \n",
    "        try:\n",
    "            ds.add_stationxml(sta_inv) \n",
    "        except Exception: \n",
    "            pass   \n",
    "\n",
    "        try:\n",
    "            # get data\n",
    "            tr = client.get_waveforms(network=net[ista],station=sta[ista],\\\n",
    "                channel=chan[ista],location='*',starttime=s1,endtime=s2)\n",
    "        except Exception as e:\n",
    "            print(e,'for',sta[ista]);continue\n",
    "            \n",
    "        # preprocess to clean data  \n",
    "        print('working on station '+sta[ista])\n",
    "        tr = noise_module.preprocess_raw(tr,sta_inv,prepro_para,date_info)\n",
    "\n",
    "        if len(tr):\n",
    "            new_tags = '{0:s}_00'.format(chan[ista].lower())\n",
    "            ds.add_waveforms(tr,tag=new_tags)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.Setup basic parameters for data-processing and cross correlation\n",
    "\n",
    "Next we setup the parameters used for processing the noise data and for later cross correlation computation. As you can find from section below, there are many parameters needed for the computation, which are associated with the input data, options for different processing procedures, cross correlation methods and some tuning parameters. Brief descriptions of each parameter are followed by the variable assignment, but note that more details on this can be found in documentations.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfile = glob.glob('./*h5')      # find sac files\n",
    "if not len(sfile):\n",
    "    raise ValueError('Abort! At least 2 sac files are needed!')\n",
    "outpath = './'                     # output dir\n",
    "\n",
    "# parameters of fft_cc\n",
    "cc_len    = 1800                    # window length (sec) to cut daily data into small segments\n",
    "step      = 450                     # overlapping (sec) between the sliding window\n",
    "smooth_N  = 10                      # number of points to be smoothed for running-mean average (time-domain)\n",
    "dt        = 1/samp_freq             # sampling time intervals of the data: in real case it reads from data directly\n",
    "inc_hours = 24                      # basic length (hour) of the continous noise data        \n",
    "\n",
    "freq_norm   = 'rma'                 # rma-> running mean average for frequency-domain normalization\n",
    "time_norm   = 'no'                  # no-> no time-domain normalization; other options are 'rma' for running-mean and 'one-bit'\n",
    "cc_method   = 'xcorr'               # xcorr-> pure cross correlation; other option is 'decon'\n",
    "substack       = False              # sub-stack daily cross-correlation or not\n",
    "substack_len   = cc_len             # how long to stack over: need to be multiples of cc_len\n",
    "smoothspect_N  = 10                 # number of points to be smoothed for running-mean average (freq-domain)\n",
    "\n",
    "# cross-correlation parameters\n",
    "maxlag       = 200                  # time lag (sec) for the cross correlation functions\n",
    "max_over_std = 10                   # amplitude therahold to remove segments of spurious phases \n",
    "\n",
    "# group parameters into a dict\n",
    "fc_para={'samp_freq':samp_freq,'dt':dt,'cc_len':cc_len,'step':step,'freq_norm':freq_norm,'time_norm':time_norm,\\\n",
    "    'cc_method':cc_method,'maxlag':maxlag,'max_over_std':max_over_std,'inc_hours':inc_hours,'smooth_N':smooth_N,\\\n",
    "    'freqmin':freqmin,'freqmax':freqmax,'smoothspect_N':smoothspect_N,'substack':substack,\\\n",
    "    'substack_len':substack_len}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Load source data\n",
    "\n",
    "We take the first station in the ASDF file as the source and the other as the receiver. The raw noise data is stored in the $waveform$ structure of the ASDF file and it uses $network.station$ and $channel\\_location$ as the two tags for each station. You would notice that most of the following steps are similar to the ones in the other notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "source of TA.K62A and receiver TA.K63A\n"
     ]
    }
   ],
   "source": [
    "# read source and some meta info\n",
    "with pyasdf.ASDFDataSet(sfile[0],mode='r') as ds:\n",
    "    \n",
    "    # station list\n",
    "    sta_list = ds.waveforms.list()\n",
    "    print('source of %s and receiver %s'%(sta_list[0],sta_list[1]))\n",
    "    \n",
    "    # channels for each station\n",
    "    all_tags = ds.waveforms[sta_list[0]].get_waveform_tags()\n",
    "    \n",
    "    # get the source trace\n",
    "    tr_source = ds.waveforms[sta_list[0]][all_tags[0]]\n",
    "    \n",
    "    # read inventory\n",
    "    inv1 = ds.waveforms[sta_list[0]]['StationXML']\n",
    "    \n",
    "    # read station info from inventory\n",
    "    ssta,snet,slon,slat,elv,loc = noise_module.sta_info_from_inv(inv1)\n",
    "\n",
    "# cut source traces into small segments and make statistics\n",
    "trace_stdS,dataS_t,dataS = noise_module.cut_trace_make_statis(fc_para,tr_source)\n",
    "\n",
    "# do fft to freq-domain\n",
    "source_white = noise_module.noise_processing(fc_para,dataS)\n",
    "source_white = np.conjugate(source_white)\n",
    "\n",
    "# num of frequency data\n",
    "Nfft = source_white.shape[1];Nfft2 = Nfft//2\n",
    "\n",
    "# find the right index of good signals\n",
    "sou_ind = np.where((trace_stdS<max_over_std)&(trace_stdS>0)&(np.isnan(trace_stdS)==0))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Load receiver data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read receiver and some meta info\n",
    "with pyasdf.ASDFDataSet(sfile[0],mode='r') as ds:\n",
    "    \n",
    "    # channels for each station\n",
    "    all_tags = ds.waveforms[sta_list[1]].get_waveform_tags()\n",
    "    \n",
    "    # get the source trace\n",
    "    tr_receiver = ds.waveforms[sta_list[1]][all_tags[0]]\n",
    "    \n",
    "    # read inventory\n",
    "    inv1 = ds.waveforms[sta_list[1]]['StationXML']\n",
    "    \n",
    "    # read station info from inventory\n",
    "    rsta,rnet,rlon,rlat,elv,loc = noise_module.sta_info_from_inv(inv1)\n",
    "\n",
    "# work out distance between source and receiver\n",
    "dist,azi,baz = obspy.geodetics.base.gps2dist_azimuth(slat,slon,rlat,rlon)\n",
    "\n",
    "# cut source traces into small segments and make statistics\n",
    "trace_stdR,dataR_t,dataR = noise_module.cut_trace_make_statis(fc_para,tr_receiver)\n",
    "\n",
    "# do fft to freq-domain\n",
    "receiver_white = noise_module.noise_processing(fc_para,dataR)\n",
    "\n",
    "# find the right index of good signals\n",
    "rec_ind = np.where((trace_stdR<max_over_std)&(trace_stdR>0)&(np.isnan(trace_stdR)==0))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Perform cross correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd5hdVbm43++06ZOZJJNOSCGUhM4AgvQSQMTQRBExXuGC7Xf1WhAUpQrYFZSrXMCLeL1UxSBoCCX0lgQSCBDSe5lM2vR21u+PtfY++5zZ58wJM5N2vvd55pmz915777Xb+tZX1rfEGIOiKIpSuER2dgUURVGUnYsKAkVRlAJHBYGiKEqBo4JAURSlwFFBoCiKUuCoIFAURSlwVBAUKCLyRRF5qRf7/1NEpvZlndxxS0TkcRHZKiIP9/Xxezj3fBE5qR+OO1NELu/r4yr9S2+/kd0JFQRKj4jI9SLy5+A6Y8xZxpj7+uF0FwJDgUHGmE/3w/EBEJH/EZGbg+uMMZOMMTP765wfhd1BiIjIMhE5LbD8WRHZLCInuuWoiNwsImtEpEFE3hKRKrdtqojMFpFtIrJKRH4qIrGQc8x0xyzacVeWHREZIyLGq6tY7hCRD0RkpFtXIyJ/cZ2azSLyv4H9fyoiK911LxeR74eco1xEGkXkn/19PSoIshD2Mu6KuBcw0tO63Yi9gQ+NMZ07uyLK9uO0xN8BZxtjnnerbwCOBY4BKoFLgVa3rRT4JjAYOBo4FfhOxjHHAMcDBvhUv17AR8B9a38ATgJONMasdpv+CqwDRgNDgJ8HdrsH2N8YU4m9N5eIyPkZh74AaANOF5Fh/XcFgDGmoP6AvdwDqgPqgd+69V8EXgZ+5dbfjBWU1wLLgQ3An4ABrnwx8GdXdgvwJjA0cKwlQAOwFLgkS12iwPeBxa7sbGAvt+1Yd8yt7v+xgf1mAj929W0B9smybgD2hVsLrHbXFA3U8aXAMX8DrAS2uXoc79afCbQDHUAjMDdQh8vd71z3aQz2A54KrAA2Aj/Icj9uyDjXZcD1wJ8DZbzjxQL1uMlddwPwFDA4UP444BX3jFa6677CnaPdnedxV3YZcJr7XQT8Gljj/n4NFLltJwGrgG+7610L/FuOd24mcCvwhru/fwcGBrZ/LFDHucBJbv2PgS5so9kI/Nbdozvc9jjQBPzMLZe4sgNzHddt6/HdwDZcm7Hv8Fk5rm8ZcBpwpXu+tYFt1a7u4/P8Pr/lPY/Auh+55/tL4B897P9vwPvuXVgCXBnYlvO5AYOAae4ZveHeq5eynGcM9j0sAu4D5mC1WG/7ZHdfonlc80jgHeCqjPXPundgDvCdfm0X+/Pgu9oftuGdi23sy7CN+XFu2xeBTuD/ATH3UX0JWASMA8qxAuR+V/5K4HFsjyYKHIHt7ZS5F2k/V244MClLfb7rXoD9AAEOcS/jQPcBXurqcrFbHuT2m4ltVCe57fEs6/6G7amUYXskb3gfBt0FwefduWPuQ1kHFLtt1xNojAN18ARBrvvkfTD/7e7pIdhezgFZ7knauUKWveMFBcFiYF93/JnAbW7b3tgG4WJ3PwYBh7pt/wPcnHHuZaQEwY3Aa+6+1WAb1JvctpOw78qN7rifAJqB6izXNBPb2B7onsWj3jVhG4F6d4wIcLpbrsm8z275FOAd9/tYd+2vB7bNzfO4Pb0bHcC/Y9/tr2CFoWS5vmXumtYDh2RsOwEriL6Hfac+BL6W4xt9zHt+gXWLgK9iv7EOXIcry/5nA+Ox39OJ7rkcns9zAx4AHnL35ED3zHoSBI+496QqY/uPgOmkOotvYrWFYJmrsULSYIXWqMC2vYEkMBH7Pc7r17axPw++q/1hVdM6XCOSse2LwIqMdc8AXw0s7+dexBi28XsFODhjnzL34l8AlPRQnwXAlJD1lwJvZKx7Ffii+z0TuDFje9o6rJ29LVgHbIP4XOB6Q19yt32z91HTsyDIdZ+8Dyb4kr8BfDbLedPOFbLsHS8oCK4NbP8q8C/3+xrgb1nO8z/kFgSLgU8Etp0BLHO/T8JqXbHA9g3Ax7KcayaBxg37cbdjG9nv4YRmYPt0YGrmfXbLXq9/ELYh+T62l1uO1RZud+WyHjfPd2NRYFupu+fDslzfMlKaTiRj2+fcvve4uh+M/QZPDznOl9y1ZGp0Hd464APgP7fjm38M+EZPz809iw6sucbbdgs9C4JtwLdDtt/ltl+GFTqfxbYLgzPKCXCYe3YVgfXXAm+73yOxmuFh+V739v7trnbkj8pewHKT3f68MmN5BNbc4bEc27gNBe7HflgPOCfYT0UkboxpAj4DfBlYKyJPiMj+OeqzOGR95nm9c4/MUdfMdXtjX8C1IrJFRLZge4BDwioiIt8RkfedY2sL1nQwOEu9e6pv8D55rAv8bsY2XH1FtmNnu7/5EHZNIwLL9RnvUU/XFHw2y7HPZjD2OX3ae0bu3h+H1SS7YYxpAWZhe7snAM9jOyQfd+s8u3yu4+bzbvj31BjT7H7mur6vYLWyu0VEAutb3P8bjTEtxph52J73J4I7i8i5WPPZWcaYjYFNU4GnAuv+4taFIiJnichrIrLJXdcnSH+Psz23Guw7m/mceuKTwHUi8qWM9S3YjsM9xpgOY8wD7tgfDxYylrdc+RsCm74A/K8rsxr7XLNed28pNEGwEhidwxFsMpbXYD8aj9FY1XK9e7g3GGMmYlX0T2IfHsaY6caY07Ef3QdYs0i2+owPWZ95Xu/cqwPLmXXNXLcS2+sbbIypcn+VxphJmTuJyPHAVcBFWDW5Cuub8D7osHPlqq9/n3rYLx+asD1Sj+1xmmW7v/DRrmnNdpw7k70yjtWBtaevxPbcqwJ/ZcaY23LU83msGegwrMnheazGchTwgiuT67h5vxvbwXqso/d44M7A+nkh15F2TSJyJvYbOccY805gfQn2nTxRRNaJyDrgP4FDROSQzAq4iKJHsb6Noe49fpLUe5yLOuw7m/mceuIV4BzgNyLyucD6eXR/drneuRjuXRWRY4EJwDWB6z4a+Fx/BbEUmiB4A+sguk1EykSkWEQ+nqP8/wH/KSJjRaQcqyo+aIzpFJGTReQgEYli1cMOICkiQ0VkioiUYT+2RqytL4y7gZtEZIKL9DlYRAZhX959ReRzIhITkc9gzQn/yPdCjTFrsY7TX4hIpYhERGS8F9KXQQX2I6gDYiLyI6y/w2M9MCZHJFLW+5RvfXPwNnCCiIwWkQFYc0++/C9wmohc5O7jIBE51G1bj/VpZOP/gGtdCOBgrM33zznK98TnRWSiiJRibdSPGGO63DHPEZEzxIZZFovISSIyKkc9n8d2Ot4zxrTjzEfAUmNMnSuT9bjb+W7kjTFmDVYYnCkiv3LrFgMvAj8QkSIROQBrJvkHgIicgn1OFxhj3sg45LlYk8hE4FD3d4A73hdCqpDAOm/rgE4ROQvrtM2n7l1Y39b1IlIqIhPJswdubHTU+cBdInKBW/03oFpseGxURC4ERgEvu/t9pYhUu+/+KOBrWBMr7rwzMq77QKxp7ax86rS9FJQgcA/7HGxEzQqsPfIzOXa5F2sCegEbOdGKdSaD7Zk+ghUC72M/zvux9/Rb2N7jJqy6/pUsx/8l1jn1lDvOPVi7bT1Ww/g21tF0FfDJDJU5H76A/Tjew9r8HyHc5DAd+BfWkbfcXWdQRfYGdtWLyJyQ/XPdp15hjJkBPIjtYc1m+4ThCqxp4NvYZ/E21lkN9l5PdKaRx0J2vxlrgpmHdejPces+Kvdj/RLrsEEK/+HquBKYgrX112Hv+3dJfZu/AS4UG4d+u1v3CrZR8Hr/72Hvubecz3HzfTe2C3fPT3F1vtWtvhirXdUDTwA/NMZ4jd4PsWbIJ8XGzAfj5qcCfzTGrDDGrPP+sNFTl2T2jo0xDdj7+pC7ps9ho4Dy5etYM9E67LP643Zc9wxsW3KfiJxjjNmEDXX9Dla7vhrrD/S+4fNIRQv+GbgDuENEirFa0B3BazbGLMW+Q/1iHhLnjFAURVEKlILSCBRFUZTuqCBQFEUpcFQQKIqiFDgqCBRFUQqc3SKxWiaDBw82Y8aM2dnVUBRF2a2YPXv2RmNMTeb6PhEEbkDIb7DDtO8ODIbxthdhE5EdgQ0h+4wxZpnbdjB2VGMlNt7+SGNMKzkYM2YMs2bN6ouqK4qiFAwiEjpautemITeg6nfYgQ4TgYvdYIwglwGbjTH7YBO+/cTtG8PG0H7ZjWo8CTswS1EURdlB9IWP4ChsgqolbpTjA9iBLEGmYFO1gh24cqrLRzIZm1VvLoAxpt4N+lIURVF2EH0hCEaSPgp1FenJ0dLKuLQDW7HZE/cFjIhMF5E5InJVtpOIyBUiMktEZtXV1WUrpiiKomwnOztqKIbNiHiJ+3+eiJwaVtAYc5cxptYYU1tT083XoSiKonxE+kIQrCY9Y98o0rNkppVxfoEBWKfxKuAFY8xGl+r2SeDwPqiToiiKkid9IQjeBCa4zJMJbGbBzERP00glS7oQeNbYJEfTgYNctr8YNkHbe31QJ0VRFCVPeh0+6lIyfx3bqEeBe40x80XkRmCWMWYaNtPj/SKyCJsF8rNu380i8kusMDHAk8aYJ3pbJ0VRFCV/dsvso7W1tUbHESg7ipWbmpm/ZhtnHrg9c+Ioyq6HiMw2xtRmrt/ZzmJF2eX59kNz+fKfZ7N0Y9POroqi9AsqCBQlT/706rKdXQVF6RdUEChKniza0Lizq6Ao/YIKAkXpgbrGNgC6krufP01R8kEFgaL0QF2DCgJlz0YFgaLkoLm9k8a2TgCSu2GEnaLkgwoCRcnBxoZ2/3enagTKHooKAkXJQV2jnRojIpBUQaDsoeyWM5QpSn+ztaWDLc3tbGqy02MMLi+iS01Dyh6KagSKEsKnfvsSJ/5sJh1dSQBKE1HcT0XZ41BBoCghLK9vBvAFQVEsSldSJYGyZ6KCQFFy0NFlzUHF8YiGjyp7LCoIFCUHnZ5GEI+ickDZU1FBoCg5aPdNQ6oRKHsuKggUJQct7V0AFMejKgiUPRYVBIqSgyYnCFQjUPZk+kQQiMiZIrJARBaJyNUh24tE5EG3/XURGZOxfbSINIrId/qiPorSV7S02/QSxfGojiNQ9lh6LQhEJAr8DjgLmAhcLCITM4pdBmw2xuwD/Ar4Scb2XwL/7G1dFKWvaVaNQCkA+kIjOApYZIxZYoxpBx4ApmSUmQLc534/ApwqIgIgIucCS4H5fVAXRelTWtq7iEWEWERUECh7LH0hCEYCKwPLq9y60DLGmE5gKzBIRMqB7wE39HQSEblCRGaJyKy6uro+qLai9ExzexexqBCNRDTXkLLHsrOdxdcDvzLG9Dj1kzHmLmNMrTGmtqampv9rpihAc0cX8UiEaAT1ESh7LH0hCFYDewWWR7l1oWVEJAYMAOqBo4Gfisgy4JvA90Xk631QJ0XJybKNTYy5+gnmrtwSuj0WEcA6i+OxCJGI5J2GetnGJr75wFu0d2pKCmX3oC8EwZvABBEZKyIJ4LPAtIwy04Cp7veFwLPGcrwxZowxZgzwa+AWY8xv+6BOipKTaXPXAPDku2tDt0edIGh2PoKoSN6moasencdjb69h1vJNfVNZRelnep2G2hjT6Xrx04EocK8xZr6I3AjMMsZMA+4B7heRRcAmrLBQlJ3G2q12noGhFcWh2+PRCG2dSVrau4hHI9ZZnKdpyNMEimI72/KqKPnRJ/MRGGOeBJ7MWPejwO9W4NM9HOP6vqiLouTD2q0tQPZ5iIMaQXHcmoaMsZPTRNy2bHiCIB5VQaDsHuibqhQka7dYjcCbjziTmC8IOolFI0RttHNeWoGXn0jDTZXdBZ2hTClIVm628w00ZRMEUecs7rCmIU8L6Eoa4tHwY/5yxofMX73Vn8PAS2GtKLs6KgiUgiOZNP6I4ab2bBqBVZY7ugzxqPgaQjKHRnD7MwsBGD7A+h06dUozZTdBTUNKwdEWCOtsbOsKLRMN+AHi0Yi/nI+5x9MI2lUQKLsJKgiUgqM5oAVkNQ0FBEEsIkQkf0HgCRo1DSm7CyoIlIKjpSOlBWRzFveFRtChGoGym6CCQCk4WoOCoDUfQSApQZBP1FCnCgJl90IFgVJwtLTbBrosEc3qLA6OAYhtp0bgFdEUE8ruggoCpeDwfAQ1FUVZfQRBjSARHEewHWMD1Eeg7C6oIFAKDs9HMLi8qMcBZWDHFHjjCJLb0clX05Cyu6CCQCk4PB9BTUURrR3J0Hh/b0AZ2DEFse3wEXioIFB2F1QQKAWHN5hscHkRkJqgPogQMA3FJDCyOP/GXccRKLsLKgiUgqMloBFAeAipIdXzj0WCPoL8z9PRmZ/2YHTCG2Uno4JAKThaMjWCEEEQ9AnbcQT29/Y5i3uWGhf9/lU+duszeR9TUfoDzTWkFBytvrM4AWQZVJYmCOycxZA711Am+QiCN5bZyWuMMYjkTm+tKP2FagRKweHNOlZVagVBuEYQMA1FxdcI8p2uErbPR7DGTZSjKDuDPhEEInKmiCwQkUUicnXI9iIRedBtf11Exrj1p4vIbBF5x/0/pS/qoyi5aOnooiQRpbzIKsRhgiDY3Mejke3KNeTRuR3jCBas25Z3WUXpa3otCEQkCvwOOAuYCFwsIhMzil0GbDbG7AP8CviJW78ROMcYcxB2TuP7e1sfRemJ1o4uSuIpQRCWgTSoEdipKvvHNJRw01l+sK4h7+MqSl/TFxrBUcAiY8wSY0w78AAwJaPMFOA+9/sR4FQREWPMW8aYNW79fKBERIr6oE6KkpXmdqsRlBXZGWZCNYJAex+LCJGP4CzOxzTkRSNtbGjP+7iK0tf0hSAYCawMLK9y60LLGGM6ga3AoIwyFwBzjDFtYScRkStEZJaIzKqrq+uDaiuFSku71QjKfI0gt2lozOCyfkkx0ZU0fiirDj5Tdia7hLNYRCZhzUVXZitjjLnLGFNrjKmtqanZcZVT9jg8H0FRzI4YDhUExnDSfjUsu+1szpg0bLuSznl09JB0LpjwTgWBsjPpC0GwGtgrsDzKrQstIyIxYABQ75ZHAX8DvmCMWdwH9VGUnHgagYhQVhTLahqKBMI5tycNNUBJPNpj494c8E1ky1S6clMz33zgLba2dOR1XkX5KPSFIHgTmCAiY0UkAXwWmJZRZhrWGQxwIfCsMcaISBXwBHC1MeblPqiLovRIi3MWA5QXxUI1gqQxBKP6fUGQZyRQWVG0Rx9B8LzZyj7xzloee3sNtzzxfl7nVZSPQq8FgbP5fx2YDrwPPGSMmS8iN4rIp1yxe4BBIrII+BbghZh+HdgH+JGIvO3+hvS2ToqSC880BFYQZNMIggO8/PDRPDWC0kSsR40geN5sZRNuAMMT76zN67yK8lHok5HFxpgngScz1v0o8LsV+HTIfjcDN/dFHRQlX1rbUxpBWVGUpizho8GBvlE/DXW4IMjMF1SaiPboLE73EYSX9bQG9SEo/cku4SxWlEwee2s1X/zjG/1y7OaARlBWFKMhy5wEgSkJekxDnelEtoKgJ43ACqBENJLVR9DQ2hF6fEXpSzTXkLJL8s0H3wb6JwdPS1AjSMRYF5LewfoIAqahHqKGMlNPlBXFqGsIjYT28UxDA0rjWX0EnkbQmTR0diWJRbXvpvQ9+lYpuzR9ndM/mTS0dSZ9jSAWldD8QdZHkFruaRzBR9IInGloYGkia9ltrSltpTWgNbR3JrnvlWVqMlL6BBUEyi5NS8ikMR5rt7Zw58xF25XP3xvA5WkE8WiEzpDJZpLGhIeP5qkRFMfz8BG43n5VaTxrg94YFAQdqXvxl9eXc920+fzp1eU5z6Eo+aCCQNmlac4hCL7xwNv89F8LmLtqK1N++xLvrt7a4/F8QeBpBBEJTQ5nAMKcxXn6CHLZ/f26tNvtlSXxrGWDIaZBQeDNqrZqc3POcyhKPqggUHZpmtvDHbkAm5tsfp7nPtjA3FVb+cm/PujxeJ6G4WkEsWgkvOeeZUBZtjTUmVpFWVEsZ90B2jptOuxc2oPnLAZo7UidY5tbv6VZB5r1REt7V/icE4qPCgJllyaXRuC102u3tgCpxjoXmRpBPCpZTUPBo3lCIVv4aKZGUFlso5GylQdr50/EIsSjkl0jaO1kQEkcSNcI1myxDu7FdY1Zj69YTvjZcxx43fSdXY1dGhUEyi5NTkHgmurl9dY84qWKzkU3jSASyWoaioSYhrL6CDKOUVkSxxhozKEVtHUmKYpFKIpFsvoIGlo7/ZnU2jqDgsAKv4XrG3f4nMfrt7WypXn3yZbaU/SWooJA2cXJ5Sz27PULN9hecTyaarnbOrt801Ha8bo5iyW0EU5mhK2mcg2F16W7RmB78VubO/xGO5OURhAJjY5KJg2N7Z3+3MpB05B3zJaOLv+adhRH3/IMx/3kuR16zp3JG0s38cKHe3bGYxUEyi5HsIcb1Ajuf3UZC9enJnDZ4hKxbXINftCm/29/fJPDbprR7djdnMX5ho/6GkF4zz3zGJUldojOnTMXc+xtz/L+2u4zkLV3JSmKRYlHI6GZSps7ujAGaio8QWDr3tzeyfptrQwqc3Mut+54+/fuYnPv7IPw2ov+8CpfuLd/Bjdmsq21g98+u3CHhwWrICgQFqxr4A/PL2blpl0/yiTY+HsO19aOLn749/l85q7XANtb3pTR4w+mbHhlcX3osX3TUCJlGupKmm7mFWNIG1CWGkcQXudsGsEjs+1UHeu2dR+01tbZ5WsEYc5iz1HsaQRtTlj86911JA2cddAwIH2sQX/TFw3rjqQ+RCvclbnr+SX8/KkPeXT2qh16XhUE/UB7Z5JfP/1hv/Sa6hrauOOZhTmdkA/NWtmtwf/5Uwu49Z8f8Pvn88v03dTWyYKdNH1iMBmb14P3TCHets3N7d0a34aQBjGzZ5XpI/DMSZkNsTEmi48gm0aQvr7SOXi94yZCRgS3Ox9BImZNQ5nCyOvpez4CTyN4dM4qRg8s5eT9bH7GHdk7zxS+O5uupOE7D8/NGjocHDXeUzhvGNn8L/96dx1jrn4i1PwI8NaKzXmFM2fiTV26ZGPTdu/bG1QQ9AMPz17Jr59eyB/ybHS3h2v+Oo9fzPiQt1ZuBqz98u4Xl/jbO7qSXPXIPM67Mz2rd31jm/uf34d8xf2zOOPXL6Q5KLeX/3l5Kfe/uix0W31jW1oUTJBgw+ZpB6udIPDmGd7gHIBeRA3AtpCc/Zn280zTUNQ5mDMb8mSGaSgeFUSyNybZNILM6wjS5nwEiSzCyMuBlDINJVmzpYVXFtdz/uEjqXDnyMc09Le3VvHG0k09lgtj0YYGXl60EYC6xpTjtb/yHy1Y18Bvnl6YV9kldY08MnsV33jgrdDt6wOaWC5/UzaC2lZQKNw5c5Gt6/rwztJ5d77CJ+94abvPVxy37+PqLH6l/kIFQT+w3vVC+jZDjmXlJvuCrN/WRkdXkov+8Co3P/G+/1F6k51szGjwvXjzTXlGe7y8yJpW1m5JN2m8vqSeyb96vscYeYDrH3+PH/59frf1rR1dHHHz01z1yLzQ/YLZQL0G1NMIyoutIFjhNJ4rThjnl/U0gqC2lPnxZ9MIMm38hvRcQyJCUSySluYhSDYfQeo6ut+vts4kiag1DUF37aXB1whSPoLH567BGLjg8FG+UGxsyz2WYEV9M//54Ny0JH6Pz12TdzTNab98gUvufh1If68yI4faO5N9Yts+49cv8KunP0wbQxHklzM+5B/z1pBMGhbX2Z5zPEPjuviu1zj79hdZH7jG5o7t15yCgqShrZNFLjDB+97Wh5j8wrj2sXd4cWHPDudG9+6vCmj09760NK2z1x8UpCBYuL6Bv721/Ta4rS0dfPuhub7JpL6xjWlz13RTH73ZpEoS2XP63TlzEVf8aRZL6hr5zdML+cYDb7Fua2uPNnzvZV5W38TzC1IvlvdRN2VpoD0BsL1hfyszRq7+dPoCPlzfyNsrtgCwvL6J15aE2+M9Wju60nqP/3p3HQDPLdiQVu6JeWu59rF3aAg0bC3uelY7geQ14Muc6vy5o0b7Zb1BVkG7cDdB0JEZPuoEQTfTEGRGoxbHo1m1mMzesddIZ6sHuPDReNQ3B2Q2ol5P39cIOrtYVt/M4PIi9hpYSoUTipk+giV1jfzwsXf9Ot3/2jIAhlQU0ZU0rKhv5v/931t8+c+zQ68lG60dXWwMNKxBM1EyaTji5hlMzeJUnXrvG1w/bT6vLNrIfa8sy3qOD9alnOqNbZ3MeG99ml+iK2m4/ZmFfP0vb3H94/P9ayh1Gp7Hq0vqmb9mG1sD73uuUGSAPzy/mAffXJG2LtjQX/PoO5z2y+dZv63Vf19WbU713I0xPPDGim7CYUtzO39+bQWX3hN+b95euYXTfvk821o7/Gf+4fpGxlz9BPe/uowb//EeN/fzxEQFmX108q9fwBg49YCh3VT4XHz34bk89d56KktiXH3W/hxx89MAjKwq5oi9B/rlNroPxLNnb23uYMH6Bo4amypz94tL2dTUziF7VXHPS0vY1trJ399eA8C86yeH1qu9M+nbPJdvbKY0nnr5V29pZtiA4m49z9VbWrj8vlkpjaApd+/x59MXpKVRWLmphW2tHX59hg8otus3N2OM4bRfPk9Hl2HprZ/ImiX0B397l0fnrGL2tadRXZrwBcE+Q8rTyn3tL3MAOGbcYH+dbxranAqXBFhW38zAsgTVZQne/tHp/PeLS/jdc4vp6EqmfYiZH39LRxeJaMTP4un9z3SCJjNzTADFseyCIFOQxKIRF5pqQusBLnw0oBFkhpB6PeJBgfDRzU3tVJfaZ+EJgkzT0Nf/8hbvrd3G5z+2N/sNq2CZG2dRHI/yjQfe4h/z7CQ32+sDWr2lhY2N4YLg8XlraGjtDHXSJ5OGVxfX8/yHdfyPEwJTjx0Teo4P1qbq9Mz7G7j2sXe58sRxXHPWAUBKEwTS8ixlS064OqDRbmyw5shkEgZXJBg+oCSt7K3/tCPTP3NkqnMR9DF4HZc1W1rY0GDX/2z6AoZVFnPBEaN4b+02rv7rO5wxaWiqXp3JrIELHrc88T6LNtjOlafdeQZw9rYAACAASURBVO/5L2Z86Jdr6+yiKBYNPUZvKTiNYGNjG147N3v55tAyKzc1+414U1sn7Z1JkknDC061W7i+0R/EBLDWvSxbWzp4c9kmX63b2tLBjPfWc8iNT3HRH1717X7GGP/467a2dhsR+8AbK9jn+0/y9sotaeuXbmzyG5al9U2+nRxsD+IXTy1Im9vWGMNfXl/uhy6WxKNsaW7P6gBLJg2/fW4Rd85M+TYemb2Sg69/yrcRezb5P726nEnXTffrE7QdZ/LoHKt9HXHz09z6z/f9jzmbv2L+mpSTrTkjp47nB1he38Teg0oBqCpNUOMay8bWzrSPtyXDHNDS3uXbYSHgLO5m7053FoO13wZj+YOE2cuDNv+wWP/2zi6K4hHfkZzpf/B8JRXFMYpiEdo6utjc3E61Cxst801DnaH7ee+C9558sK7BFwKQGoexaEMD//F/b3UzFf1j3pq0Z7FmS4v/roN12IMd2T3NdWJGVqU3rttaO3hk9irau5JpvfZMgdrRlWRTU3vad+XVO6j5ZhNeG7al6h40Kc1blfqGPnPXa5x9+0uc89uXOObWZ7n/teV89+G5QLr9/6WFG/neI/OYv2arP0YlyOK6JjYHUnt82x1j1jLbngTfhVuefJ+v/u+c0Lp5eH64ba0d3Z5lNNC5WtqPDuQ+EQQicqaILBCRRSJydcj2IhF50G1/XUTGBLZd49YvEJEz+qI+uQg6zN50v5dtbOKMX73AL55awMbGNo7/6XN8/2/vADDpuul88Y9vsHJzM60dSWIR4aVFGzn/zlf84yzb2ERTWyf/ft8sPv37V/lwvX15trZ08C2XVx9gqbNnNrV3+aGAq7e0+PHwHrf+8wM6k4Y7nkl3mC3daI+779ByVm1qZkNDG5WuV/j9v73DHc8u4v5AL2lzs/0IPcYOLqMzaXwn5OtL6jndqaTzVm3hh39/N+18VaVx5jgT0DsuAsJ7Ueev2ZbWy53+7jrOvv1FX+3P5lT97xeX8p4TTHUNbRhjWLShIc055p2rIpCvx/sItrV2YoxheX0zYwaV+fuUO42lobUzrafaTSNoT01KA6nRyGEaQaaCk8s0FJam4msnj+fasw8gGpGsPoKiaIR4LIuz2PX0yxIx/9ybm1MaQTwaoSQe7da4eNfi9Vo3ZvEFtHZ0sXRjE5fe8wbT5q7hn++uZWtzB5fc/RpL6hr5+l/e4uzbUw7Pt1ds4ZHZqzho5ADAmuDaO5Mcc+uzPPOB7S2nzcPcmeTjtz7LVY9aX9At5x3EKfvbSKeNGR2HW5/8gMNvmsHbK1OdM88vtGRjE11JQzJp+DDDOVtVGqcsEWVjY5t/3cGOwPw13cdvePzwsXd5ePYq/uflpf43C/D5e17nwVkrOfv2l7jrhe62+WAdAcYMKuWR2au4bpr1hwXHszw1f13GvumdO0iFBa/d0tot8i2onS8KEUp9Ra8FgYhEgd8BZwETgYtFZGJGscuAzcaYfYBfAT9x+07ETnY/CTgTuNMdr9/wejRjBpXy5rJNrNzUzGfveo3FdY3c8ewivuJsjm+v3OK/1K8sruc990Kdf/hIIP2F//lTH3Liz57jjWVWsHi9v+X1TTS0dXL5cWMB29M9+pan+eucVOO8YF0DmR10b3nhhkZ+Nv0Df1Tj0o22t3TE3tVsbu5gQ0MbYwaXUVUa9/cJ9lQeeHMF6wM9pXE1tuHc4sxDX/vLWyzc0MisZZv41G9f5n9fT7eP1gbMXV4D6L2ox08YnFb2jy8vY/6abVw3bT6vLq4PDWk8bp/UPgPLErR0dNHU3sW5v3uFj9/2rL9t3iorCAZXFNHc3kVDq73WgWUJupKGjY3trN7SwtjBKUHg2fxbO7vSHOJhPoKSgEktljN8NF0SFMWjWZ3FYRrBd8/Yn8uPH0dpPJrVNGQ1gqirQ3dncXlRjGhEfG1kc3MHA51GANZ5nnmvPe1mwzYraLM5hZMGPn/36zS2dlJdGueNpZt4dUk9Ly+q51sPze1W/s+vL6exrZMfn3cgYJP+1Teljr3v0HIa2zp5ZPYqxl7zBC8tqkub/W3/4RV8/mPW7JJZpz+9ugyA5xbU+Vqn5y9r70wy/vtP8ptnFrIkI7fS37/2ca75xAEkTco3tCZkoqFcXP/4e9zxbPYoJU/L8Rpsr0G+4VOTGDe4jPbOJP8MzCm9ZmuqUxM0l1aVxrnvFdtRe2PpJs6782VWbmr26716SwuNbZ1+5w7SQ6K9AI7+oC80gqOARcaYJcaYduABYEpGmSnAfe73I8CpYu/QFOABY0ybMWYpsMgdr9/YsK2VRCzC6ROHMnflVv7wwmI2N7dz35fsad906l08GvF78AB3v7QUsB/3104e768f4WzmmVE6AO864XHk2IHEo8LzH9axflsbd79ojzWsstjvCYclTFuxqZnfPbeYL9z7BsmkYenGRmoqihhZVUJLRxerNzczxC17BGPwfz9zsR+DDjDONZz/eGcNHV1Jv1cW7A15jeQJ+9bw3TP289d7vaxtLR0cM24Q9192NONrUg3xko1NFMUijB5Yytf/MidNJQe48IhR/OlLqUd76F5V/nEzGzLPpDGsspitLR0scc/B2+etFfYZ7Ts05WPwzD0t7empJcLCR4NOfM8+Hxo+SjrFsUgOjSB7KGVJIprVWZxwvgQIMw11+E7n4niUVpc2o7o09UwrimLdepHeuTY0tLG1pYP2rmSa0LzlvIP4ykn2HV69pYXTJg7l+Ak1vLF0k//+hPVc129rIyJwwPBKyoti1De1+w36F47Zm3MPG0lX0nDD4/MxhjRtFGBEVQk15fZ7qWtoo7m9k6Ubm1i3tTXt/nkjpj1H7H5DKwB4eNbKNBNkVWmc0QNLGeKc6Z55aN3W9NDLQQHBmY1cobWH7GU1IE/Ye5FKE4aWc/qkoWxsbGdxXSNnHTiM/YdV+AkBwd7fcYPLeOCKj/HFY8fw9PvreXPZJn7//GLeWrGFr/7vHP8ert3aQmNrZ9qz8u5LdWmcB99cwaIN/TO2py8EwUhgZWB5lVsXWsYY0wlsBQbluS8AInKFiMwSkVl1dR8978f6ba0MrSziyDEDae9K8pfXVzBxRCXHjh9ERSDSY+WmZhbVpW767OWb2XtQKTUVRXzr9FQDme37L45H/Bdn70Gl7FVdyuvuZfNs5JNGVPrl93UvezZ+On0BS+qaGDu4jAGuIVhc10RNRXGaIAjmtdnW2skFR4zyr2usa7h/+q8FPDV/vV/OU7fvmVrLvOsn89evHsvvPncY+w2r4MWrTmbCkHK/l7WttcMPjSyOpytvk0ZU8sd/O5L6pnZfu/BCH4dWFvnTPUKqUQ9GiQQRsftsae7wM2x6+8x2gmCfIal75msEHV2hpqH2ziRjrn6CGe+tpyTgI8geNdR9iszieJS2HqKGbpoyiRe+e3LattJELo0gSjyWzVnc6YfLlsSjNloladIEQaZGELQzr9vawg8es+Y+zzEfjQifO3o0VwbCbkdVl3DwqAFsaGjLGhfvMayymHg0QnVZnM1N7X5n4vzDR/kBBZ5m8+Q76xjj/Dhgx1YMrrB139jYzjV/fYeTfz6T3zgT6D1Ta6mpKOKSj+0N2Ea0ojjGk984nsuOG8uGhjZWbmqhypnGDh5VhYj4PpMtLfa5r9rcQkTwfUhe1FUmlxydcgpvyBFKm+n38Bru0kSMmvIi2ruSLKtvZp8h5ZQXxboNujv/8JF8bNwg/v34cYysKuE7D8/lVedAficw6Gyt6xTtHTB5evzwkxNJGpge+G77kt3GWWyMucsYU2uMqa2pqfko+/Po7FU89vYahlYUc+QYa/ZIGtuAiYj/QR6yVxVtnUleW7yJiMDPP30IFcUxfnrBwYD9mG49/yD+cOkR/suXyfiaVG91VHUpew0s7WY+CAqC/YfZRu2QUbb3EWyDzj98JL9/fjGzlm9mwpDytEFUQyqKGFkdEASuwZ5y6AgALqrdiy8caz+sw0dX++UempWSv54gqCqNE49GOHx0tT9Yaa+BpYwdXMZLC+t4d/VWGlo7/Q9+r2r7oZ20n30e1aUJxteUM6gswWKnPntRRsUu2uGYcYMAOGy0bdTfWZU++tK77pJ4lKrSBFua23l75RYS0Yhvm569bDPxqPgfOlizDdge/+bmdl9T83rHs5alenxBH0G2GP7MXEOAHUeQxVns9dyOGT+Y0YF62fPFsgwosxFMnrM4M99QY1unHxk0tLKY911ETXWghzugJJ42J0HQmfr0+xt4wjmHJzhB4AnM4Ds0qrqEEa6xC9MEguW9d21gWVGaRjC4POHXNXiPDh5VxRUnjOPEfe07MqjMNsorNjX7UXL/94btNBw2upo3f3AaU4/Z29+/ujRBNCIcPXYgnUnDik3NTBxuvxvvW6lydfPuw/trGxhXU+4HEHidkSD3X3YUXz15n9BrBTh2/CAuP24sV54wjq9lKVeWiKYJmfE15b7gDjLaNexlRTG+d9b+LK9vpqWjy/eXgP3+l25sYltLB4PKE2mdUoAJQyqYNKKS5/sp+V1fhI+uBvYKLI9y68LKrBKRGDAAqM9z3z5BRLj3ZWuSGTqgmOqyBPsPq+CDdQ1+o+05myZPHMrclVt4bsEGRlWXcuERozjvsJFp5puLXfx6toZhnyHlzF+zjerSOOVFMSYMKef5D+sQsZEAnUnDUWMHYa1hKY3g4FFVzF21FWPgpnMPZFtLB5cdN5a/zrG3ZeqxY9KcYUMqi/wPMMgPzj6A8w8fxfiacr59+n5ceeJ4KovjvHvDGRx0/XSe/9DaYkviUd80NKAkXIUuTURJGvjkHS9Rloj6QuLW8w/iiL2rOXR0FTMX1PkN7NDKYt8hPHpgKe+s3krUmT9+//kjeG1pPQeP8sw8tuG55byDeHTOKkriUV5atNEJgjjbWjv525zVnD5pKIOcmWve6q2MGVSWNojIMw15dvSR1SWs2drqm4aCYxZK4qn7Fcs6oCw915A9hzXPhOGlnggz8ZUmot2ilzq7kiQNfooJCHcWe892RFWJ3wh4zmKw93rh+o3+8mNvraY4HuGIvavTbMr7uY6G94yC2s5e1aUUu/Vvr9ySFvbqcZh7xl5dB5bGeW5BHS8utOceXF7UbewE2F730U74g02hUFUa5wEXrz+yqoTVW1ooiUf964pFI75PxOv9H+K0QYDavavZZ0g5Fxw+CkgJqa0tHTwxby3PLdjA2QcN9x30QfOoiBXyVSUJ36QUxq8+cyhDK4v95TCtriQR9YUNOEEQcg+CJsyzDxrOk/PWUjummurSBM86J/s5h4zgZ9MXANbcN6g8keZfqSqNc8K+Nfz3C0toaO3wv8G+oi80gjeBCSIyVkQSWOfvtIwy04Cp7veFwLPGxmtNAz7roorGAhOAfkvzd/pEG9/rhWR9/xM2NtlzYl55orWbTnblNjS0+T31bJOe/Pmyozn7oOHd1nvq5EGuwTvU9YCNgSPHDKQsEeXY8akPxOvdHD0u5aC99GN787WT96E4HuXuL9Ryw6cmse/QirTe3MiqEkZVp6uuiViEIRXFfi8sEhG/F19eFGOYe8EPHFnJ0Moi3zYdPG6QSSMG+L+b2rt801B1WYJ/P2EctXtXc/05E/nROTZGYNiA1Af0tZP34coTxzH1mDH2HKVxzpg0jAElcSqKY8xxZp5TDxjCo185liGV9sMqSUT9nl5DWycXBkwP7Z1JXyh4eD1dL/30kMpiP1pna0uHL0i9Y3t4z7W7RpAtfDT3OIJYFkGQ2Yh4ZiAv6VxYHewHb+/1yKrUPa0KmIZGVJWwvqGVji47qvcf89ZyxqRhHD8hpTXfdekRvvZZEu8eizGqutTX3No7k762bM9r3y1PcHsdn4Fl6Y1ocTya1jide+gIfnrBwWlCwGPMoDK2NHcwuLyIC48Y5Y6XSBNO5UXxtGsNNtrDBpRw45QDGeNs6V5epyV1TXztL3PoShoOGF7prw/Wy3uHBpTEu41G9jROIK2BD96HIKWJGIMD9Zo4ojK0gd4vYPaNRoTfX3oElx8/Lk3QeEEoYM19n67dK23cUVVpnAsOH8lvP3dYt3r3Bb3WCIwxnSLydWA6EAXuNcbMF5EbgVnGmGnAPcD9IrII2IQVFrhyDwHvAZ3A14wx/ZZc/eT9hvDrpxf6N/KEfWtYfMsn/MbgqyeN54oTxhGLiP/xHjC8MtchOW7CYI6bMJiJzy0iGhFuc4NSzj1sJA2tnXz3TOtPODTQo/n25H1ZsL6BSET4j1P24an31nPM+EH84/8d52sGYzLMC6dNTA1SqSpNV+u9uW89yhK5A6/GDi5j7dZWPnPkaB6fuwbrsskuCP7t42MoK4r5IbWZg91EhC9+fKy/HHzBh1QW+YOBMhlVXcr7a7cRkZT67qnEnmnI44DhlX5vNKyunr+ipd1GDQ0qS1DionV+/fSHbGput7H4nck0H4H3LmSa7bKHj+aOGgrrMJTEo92iZNrccYpiKdNQppBpbu+iNJHSCDwmBHqYIwYUY4z1fa3Y1MzWlg7OOnA4k0ZU+u/i5EnDaO3oYq+BJdwwZZK/7wHDK3l/7TaGDSgmInYinqQhbaCV18P1zDEnOzPgwLLu70qwN3z58eM4MNCwBvnqSeO54v7ZnLxfjR/Jlnn/K4pjbGxs87UEEaEkHqWlo6ubzb84HqU4HvHH+YCNUPJt+UX23bjk6NG8tGgjW1s6/Pfn6rP2pzQRZUV9M8fvW+OPjI5kPMe9B5WycEOjf4/ACvjRA0s5euxA/uPUCUQj4gvuhJtjIhaRrAMthw1IXcfwASUcP2EwizY0cvTYQRyyVxXHjh/EeXe+QiwilBfFqBhSkeYX60v6ZGSxMeZJ4MmMdT8K/G4FPp1l3x8DP+6LevTEwaMGcPvFh6WFMQY/XBHxIzjyFQQenh3R+/j2HVrBTece6G/3ehTja8qoHTOQWtfr+tbk/fjWZCssvA/nnqm1TByR/bzBRnBEVYnfqHiU5khtAdbnsW5bK4ePrvZt5yWBVAeZxKKRNHtmmCkqyLCAIAhTlT1GVZfw/tpt1FQU+c/BGyRVkogyICDwhlQU0RWIs80URl5Pt7HNagBVpQlKElEWrGvg9aWbuPio0azZ0mJNWMHw0WzOYrqHj+YaR+CZoHrSCJ58Zy0Pz1rJredbf1MiFqXMNVRho6C9QVhBQRC89uFu/dqtrTzz/gYSsQjHTxhMWVGMK04YxyGuJ18cj/LiVaekHf/Plx3F4rom/7l7DdyEoeV8Z/K+rN3ayvETavjyn2dzzPhBvHrNKQytsM82KKQ9gu9F5qjxIKdPHMp/XXI4x+4z2B8fkhm15b0PVYF3fWBZgtVbWtJMPR5VJQk/rPOeqbWctG8Nb7kBo0WxKMtuOxuAT/32JURSdf3yiakIwGyBC5ByOA8fUOJH+hXFIogID155jF/Oe98PGF7BeYeN5IR9s/szhwS+E4D7Lzs6bfmw0dW8eNXJlCSiWYVJX1FQKSZEhE8dMiKvsrV7D+Rf89dxwPDtk8A3nXugnx8n89wvXnVyj40o2NQXuQiqn6WJGKUJuPeLtby4cCN/fHlZjxlDR1SlnIOeBtLTLFdDK4sYO7iMZfVNOT9ywI8Mge6RRUG8j+vj41OC2XO2RUTSBF4kIkRIaWrZNIINDXbkeHVpnNJElLdWbKErabj4yNHc5RJ3hYWPZpplktZJkH4Op1FkRhQ1tHbwu+cWMXZwWVqMv0fQWeyNMvUanaJYxDdVNbtnsG6rjWxr7ejyr8sz3WRet2cyWrOlhffWbOPAEZW+MPVMn9kYVF7kp68Icv7hIxlSkWqkvEY0eG4vMunIMdV878z9gXRBkOu5iwhnOXOqp/kGNWZIxeqfeWDK7Drl0BHcOXNxaBTQgJI467a1MmFIuf/9xHxtL5lWrrI43q3Hb+uf3e7uPddYNL3jmIknCIYNKE7TksPIdAiHsdfA0h7L9AUFJQi2h599+mAuOGJUaChXLi792N5Zt/XVQw0zP5yy/1CGVZbwx5eXhY5pyEa+Go+I8PS3TqQrabJqDh6eX+XmgEYUhtd7PvPAYf4670NKGpPWG/SoLI7T3N7l2389ilydNrg8Q1Wlcd+UALaR8kxmYQPKuo0DMIQOKAMb/x9s6JbXN7OxsZ0bPnWg3/gEKU1EaWnvTBPQXv6ZRCzia3At7Z3MX7OVs29/iR+fdyCtHanzjKou5TO1e3HpMenvl2fGWbOllQ0Nrb5T+KNwy3kHsbm5PU0IZMNzgH7ztH197TaX9peNqtIED115TLcO1+0XH0ZdQxvHBPxo35m8HxfV7sWo6u7f0YBASKmH/2wD2t6wymLWV4YPOPMGcu0fcg+9aKew+a2DeJrNsMqe76GIcN5hIzl8dFWPZfsbFQRZqCiO+87l3QXvY8qnp+ER9tJnIxqRrE7zIEfsPZC5P5qcZtoJ4xunTmCfIeWcFtCAylyj2JU0oeYHr+eX2TOORGyaaC8evKIontYwlRZF/d5/SSI4jiCbRmC6DygLjFUICgJvxKln4smkNBGluaPLHzEN+LmbimIRXzA1tXXx+hJrqpuzfIs7p61fNCL85MKDux27rCjGgJK4TYS2rS3NSby9fC4QV98T5x46kkNGVTEuECYdi0Y4Y9JQPhESPJGLoFPUI0xzj0TEdxBn4r0P3uAvgHik+/iMq8/aPy3NeZCK4jh3XHxYWsCGx0VH7sXs5ZuZPGko33jg7ZC9LV4a6Xw7fb/6zKF5letvVBDspvzzG8d36xWLCK9dc2q3aJdceGaEXKF0H4WehABY08QXXDSRh2ca6koaG1lUFOOagInDi1jKzPcPtqH2BEFlSTwtpru8KOYLgmBPPx7SawQbPtrdR2D3b8uI9/d6+tkyQ5YkohiTGuxXWRzjAxfvn4hF/BQSLR1d/ihhLyAgLMonk+EDillc10hDW6cfddXfiEiaEPD4w6W1O+T8mXjaY1AjmDTSartBs5M1h2U/zjlZTMflRTF+d8nhaUn4wvjSx8ewraWDS47ObhnYFVFBsJuSzaQTDN3Ml9euOTUtI+fOxOvFG2N7we/ckJ6HMBj+l0lxPOKbhiqKY/6xROyANj+VQ6CHGAukmAim+U0a0z1qKJbSCIJ4gqEoyz300oV72VZrKorY5tIUeJpFaSLmJqW3gsxLJJePIBhRVcLrbk6IoXmYdfZEhlQWURyPpJmYjh0/mBevOrlP7ew9pa2vKk1w/acm5SyzK7JrfP3KTmXYgOJQM8zOwNNQurKkyva0oLAPsiQe9SdpqXTjFMA2xJGIhNqM4059+uuc1ex37b/8CW9MWK4h3zSUoREEQkHD8HwAXuqB4EjXikAKiUdnr2bGezaFgCc0cjldPUZUFdPknNE7SiPY1fj348fx6FeO7aaV9bWzNZ9gj90RFQTKLoXn0E1mSeLkOfTCNYJUIxDUCDzhEpZy2tMIvDxQf31rtZ+bvnuuofB4f0/DyGUaglR2zGDUS4UbOFVWFE2L3PJy+OQjCIJx/0PzcFLuiVSVJtIGPvYXZR/BIb47oIJA2aXwQ/560ghCBIEX1SMC5YmYPzrVM9l88mDrxDwr4MwMhgMCPP3eej+ld9iAMugeauslosuuEdj9NjV5eXlCNAKnNZy4bw1jB5f5kV/5mOyCo14L1TS0o+iPUb27AnumeFN2W4YPKKaiKObHpmcyYWgFg8sTaRk4PbwRw+VFMSIR8Z3FxqT29WLiPTIHgC2vb8ITQZnOYq9nn5lSuicfQYkvCNoRSU+L7PUwPT/CoLIE9U1tvmM5Hx+BlwjuyDHVeTnpld5z7qH5jUfaXVBBoOxSFMej3RzEQc45eDjnHDw8dDCP12P3/AcVRemCIIxYxgz1TW4iHOjuI/BCWzNHAPuCIItpyPMR1De1UxKPppkXvDEZXujpwLIEZYmYf8ySHtKFABwxupo7Lj5stwt33l3J7EzsCaggUHYrcg2196J6PHNLSiPILgni0e7HW+cijzJHn3omnsxpJ1Phoz2ZhqwgCEtV7Imd6rJE2viHfHwEkYhkDXtUlHzYMw1eSkHi9Z49/0FqlHL2fURSg+S81Bnrskx16KeCyNQIXNRQIov92DPvbGnuoDgeDR2B6wmTQWWJNI0hH9OQovQWFQTKHoOXXsELAfU1AnKnBfA0Bi9V83pPI8jQPnKZhhLRSGj+GkhpBGCdv2GRJ14k0sAMQZDN76AofYm+Zcoew2XHjeW0A4b6GSUr8tAIgts9QbDWaQTdo4YiiNAtqaAdiJb9Uwpmgy1JRCkPSUWRyvOfSNuuGoGyI1BBoOwxxKMR7p5ay5eOs1kfy/JwFgc5Zf+hbm5gG+aZ2b8XEUrjUX/wloedezj7p+QJELANuxfWGiSVHC+ephHk4yNQlN6igkDZY/FMOZNyzO0AcMHho7j27APYZ0g5A8sS1LvBXJmmIQiff7itM5k1YghSk6qAbdjDktOdOclmYB1amT7l454at67sWvQqakhEBgIPAmOAZcBFxpjNIeWmAte6xZuNMfeJSCnwMDAe6AIeN8Zc3Zv6KEqQAaVxHrziYzkn+QH4xUWH+L9jUfFnywoLUPJSSgexgiB3g+3No1CSxVn8n6fvy+XHj6WqNLHLpPtQCofedjeuBp4xxkwAnnHLaThhcR1wNHAUcJ2IVLvNPzfG7A8cBnxcRM7qZX0UJY2jxw3arom+hdSo5rBQ1dJEd9NQW0dXj3M0eBFHJYnUOILgiOBoRHwB4E0HqSg7it4KginAfe73fcC5IWXOAGYYYzY5bWEGcKYxptkY8xyAMaYdmAOM6mV9FKVXRCSgEYRstxpBF6s2N/uRPvloBJ6ZqiRu50X4zWcP5cErPxZadlB5EfsPq0iLNlKU/qS3gmCoMWat+70OCBvaOBJYGVhe5db5iEgVcA5WqwhFRK4QkVkiMquuri5bMUXp7HCvdAAAEtBJREFUFSI2BTUQOq9DWVGMpvZOzr79Jf5r5mKAtPTV2fCyYHrO3ymHjgydactj2tePY/a1p3+US1CU7aZHQSAiT4vIuyF/U4LljA3GzjM+I+34MeD/gNuNMUuylTPG3GWMqTXG1NbUqOqs9A8RET9NdZhpqCQeZVNTO1tbOnhntZ2kpK2HqCGAw9x0hD3NJ+2RCMxlrCj9TY/OYmPMadm2ich6ERlujFkrIsOBDSHFVgMnBZZHATMDy3cBC40xv86rxorSj0REcmoEpYmoP/L4w/V2lrG2jiSDynILgkPczFkL3MxkirIr0VvT0DRgqvs9Ffh7SJnpwGQRqXZO4sluHSJyMzAA+GYv66EofYJIYCL7MGdxUSoh3KrNLazfZieN78k0dPAomyt/e+fzVZQdQW+Tzt0GPCQilwHLgYsARKQW+LIx5nJjzCYRuQl40+1zo1s3CvgB8AEwx6nhvzXG3N3LOinKR0ZE/JHGoc7ijAFeR99i3Vo9RQ1VFMdZ+OOzdFyAskvSK0FgjKkHTg1ZPwu4PLB8L3BvRplVhH9rirLTiEhqdrSwAWWlWWaoWrShscdjqxBQdlX0zVSUABER3zQUNqBscHn4YK9oloRzirI7oIJAUQKkawTdtw8JzDdc7WYD+/QRo7jzksN3SP0UpT9QQaAoAUSEzqR1BkuI5bImMCfw4aPtAPnPHT2aEYFRwoqyu6GCQFEC2AFlqd+ZBDWCI8cOJB4Vxg0u30G1U5T+QaeqVJQAaSkmQiRBTUAQfPHYMZy4b41OGK/s9qhGoCgBIoIvCMJ8BMH5AYrjUQ4YnjuzqaLsDqggUJQAIrnTUCvKnoiahhQlQEQCaaizDHOpqSjyI4sUZU9ABYGiBIiI+I18No3g5e+dsgNrpCj9jwoCRQkQzDUU5iyGntNJKMruhr7RihIgqBHoYGGlUFBBoCgBRKRHH4Gi7GmoIFCUAD2FjyrKnogKAkUJENHwUaUAUUGgKAHSwkdVEigFggoCRUlDMDkmplGUPZFeCQIRGSgiM0RkoftfnaXcVFdmoYhMDdk+TUTe7U1dFKUvCPoFVCNQCoXeagRXA88YYyYAz7jlNERkIHAdcDRwFHBdUGCIyPlAz9M7KcoOIDgrmTqLlUKht4JgCnCf+30fcG5ImTOAGcaYTcaYzcAM4EwAESkHvgXc3Mt6KEqfEAl8EaoQKIVCbwXBUGPMWvd7HTA0pMxIYGVgeZVbB3AT8AuguacTicgVIjJLRGbV1dX1osqKkp2gOUhNQ0qh0GOKCRF5GhgWsukHwQVjjBGRvDNxicihwHhjzH+KyJieyhtj7gLuAqitrdWMX0q/EDQNqRhQCoUeBYEx5rRs20RkvYgMN8asFZHhwIaQYquBkwLLo4CZwDFArYgsc/UYIiIzjTEnoSg7iWDjH1GNQCkQemsamgZ4UUBTgb+HlJkOTBaRauckngxMN8b8lzFmhDFmDHAc8KEKAWVnkx41tPPqoSg7kt4KgtuA00VkIXCaW0ZEakXkbgBjzCasL+BN93ejW6couxzpUUMqCZTCoFdpqI0x9cCpIetnAZcHlu8F7s1xnGXAgb2pi6L0BaI+AqUA0ZHFihIgbeyASgKlQFBBoCgB1DSkFCIqCBQlQLDtVzGgFAoqCBQlQNBHENEcE0qBoIJAUQJEVCNQChAVBIoSIKIpJpQCRAWBogTQAWVKIaKCQFECiEYNKQWICgJFCaBRQ0ohooJAUQLoOAKlEFFBoCgB1EegFCIqCBQlgGoBSiGigkBRAqizWClEVBAoSgA1DSmFiAoCRQkQbPxVI1AKBRUEihIgfWTxTqyIouxAeiUIRGSgiMwQkYXuf3WWclNdmYUiMjWwPiEid4nIhyLygYhc0Jv6KEpvSfcR7MSKKMoOpLcawdXAM8aYCcAzbjkNERkIXAccDRwFXBcQGD8ANhhj9gUmAs/3sj6K0ivSG3+VBEph0FtBMAW4z/2+Dzg3pMwZwAxjzCZjzGZgBnCm2/Yl4FYAY0zSGLOxl/VRlF4RUY1AKUB6KwiGGmPWut/rgKEhZUYCKwPLq4CRIlLllm8SkTki8rCIhO2vKDuMYNsfi6gLTSkMenzTReRpEXk35G9KsJwxxgBmO84dA0YBrxhjDgdeBX6eox5XiMgsEZlVV1e3HadRlPwJTkajckApFGI9FTDGnJZtm4isF5Hhxpi1IjIc2BBSbDVwUmB5FDATqAeagb+69Q8Dl+Wox13AXQC1tbXbI3AUJW+CkUKqESiFQm/f9GmAFwU0Ffh7SJnpwGQRqXZO4snAdKdBPE5KSJwKvNfL+ihKrwj6CKLqJFAKhN4KgtuA00VkIXCaW0ZEakXkbgBjzCbgJuBN93ejWwfwPeB6EZkHXAp8u5f1UZReEWz7VRAohUKPpqFcGGPqsT35zPWzgMsDy/cC94aUWw6c0Js6KEpfohqBUoioEVRRsqCCQCkUVBAoSoCgRhBTQaAUCCoIFCWA+giUQkQFgaIECI4jiGrWOaVAUEGgKAHSks6pRqAUCCoIFCWAtv1KIaKCQFECiGYcVQoQFQSKEkA1AqUQUUGgKAF0ekqlEFFBoCgBVA4ohYgKAkUJoBqBUoioIFCUAJp5WilE9LVXlAAaNaQUIioIFCWAWoaUQkQFgaIEUB+BUoioIFCUACoIlEJEBYGiBNABZUoh0itBICIDRWSGiCx0/6uzlJvqyiwUkamB9ReLyDsiMk9E/iUig3tTH0XpLV7SOVUMlEKitxrB1cAzxpgJwDNuOQ0RGQhcBxwNHAVc5yayjwG/AU42xhwMzAO+3sv6KEqv8ASAmoiUQqK3gmAKcJ/7fR9wbkiZM4AZxphNxpjNwAzgTEDcX5nYblglsKaX9VGUXuEJADURKYVEbwXBUGPMWvd7HTA0pMxIYGVgeRUw0hjTAXwFeAcrACYC92Q7kYhcISKzRGRWXV1dL6utKOF4AkBUI1AKiB4FgYg8LSLvhvxNCZYzxhjA5HtiEYljBcFhwAisaeiabOWNMXcZY2qNMbU1NTX5nkZRtgtPI1AxoBQSsZ4KGGNOy7ZNRNaLyHBjzFoRGQ5sCCm2GjgpsDwKmAkc6o6/2B3rIUJ8DIqyI1EfgVKI9NY0NA3wooCmAn8PKTMdmOwcxNXAZLduNTBRRLzu/enA+72sj6L0CvURKIVIjxpBD9wGPCQilwHLgYsARKQW+LIx5nJjzCYRuQl40+1zozFmkyt3A/CCiHS4/b/Yy/ooSq9QjUApRHolCIwx9cCpIetnAZcHlu8F7g0p93vg972pg6L0JREdR6AUIDqyWFECiEYNKQWICgJFCaA+AqUQUUGgKAFSgkAlgVI4qCBQlACpAWU7tx6KsiNRQaAoAdRHoBQiKggUJYAnAKIqCJQCQgWBogTQ8FGlEFFBoCgBIjqgTClAVBAoSgAVAEohooJAUQL4KSb0y1AKCH3dFSWAoOMIlMJDBYGiBPA0ARUESiGhgkBRAmjUkFKIqCBQlBBUDiiFhAoCRQlg3GSrahpSCgkVBIoSIOkkgQoCpZDolSAQkYEiMkNEFrr/1VnK/UtEtojIPzLWjxWR10VkkYg8KCKJ3tRHUXqLJwhUDiiFRG81gquBZ4wxE4BnyD75/M+AS0PW/wT4lTFmH2AzcFkv66MovUJNQ0oh0ltBMAW4z/2+Dzg3rJAx5hmgIbhObHavU4BHetpfUXYUviBQo6lSQPT2dR9qjFnrfq8Dhm7HvoOALcaYTre8ChiZrbCIXCEis0RkVl1d3UerraL0gG8a0rghpYDocfJ6EXkaGBay6QfBBWOMERHTVxXLxBhzF3AXQG1tbb+dRylsPJNQSSK6k2uiKDuOHgWBMea0bNtEZL2IDDfGrBWR4cCG7Th3PVAlIjGnFYwCVm/H/orS5xw4spL/OGUfPnf03ju7Koqyw+itaWgaMNX9ngr8Pd8djTEGeA648KPsryj9gYjwrcn7MWxA8c6uiqLsMHorCG4DTheRhcBpbhkRqRWRu71CIvIi8DD/v727j5GrqsM4/n1EqYlgsIBYAe1LQMVGpTTYGNTEGoFGqYSQkJiAEUOMkkgIkZL9p4F/qFU0xBeCgQRMI4nyImKQF4NimhQtUHYX29JtqaFNbaFG2gZdMf35xzmDt5uZnZnu3Hune59PMtm759679+npmT1zX/YcWC5pl6QL86obgeslTZDuGdw1wzxmZtanrpeGphMR+4Hlbco3Al8vfP/pDvvvAM6fSQYzM5sZPyRnZtZw7gjMzBrOHYGZWcO5IzAzazh3BGZmDeeOwMys4RRx7I3WIOlV4G9HufspwGsDjDMoztUf5+qPc/Vntub6YEScOrXwmOwIZkLSxohYWneOqZyrP87VH+fqT9Ny+dKQmVnDuSMwM2u4JnYEd9YdoAPn6o9z9ce5+tOoXI27R2BmZkdq4hmBmZkVuCMwM2u4WdsRSForaYukUUkPSjqpsO4mSROSthbmRkDSRblsQtKqknJdLulFSYclLS2Uz5f0L0mb8uuOwrrzJI3lXLdLGviEup1y5XW11VebnKsl7S7U04puOatSR31Mk2VnbjObJG3MZXMlPSFpW/76ngpy3C1pn6TxQlnbHEpuz/U3KmlJxblqb1uSzpT0lKS/5vfjt3N5uXUWEbPyBXwBeHteXgOsycvnAC8Ac4AFwHbguPzaDiwEjs/bnFNCro8AHwL+ACwtlM8Hxjvs82dgGSDgUeDiCnPVWl9tcq4GbmhT3jZnhe2tlvqYJs9O4JQpZd8FVuXlVa33RMk5PgMsKbbtTjmAFbl9K7f3ZyrOVXvbAuYBS/LyicBL+fil1tmsPSOIiMcjzYUMsIE0JzLASuC+iJiMiJeBCdLkOOcDExGxIyL+A9yXtx10rs0RsbXX7fNc0O+OiA2R/ufvBb5cYa5a66sPnXJWZdjqo52VwD15+R5KaEdTRcTTwD96zLESuDeSDaQ5zedVmKuTytpWROyJiOfy8kFgM3A6JdfZrO0IpvgaqdeEVKmvFNbtymWdyqu0QNLzkv4oqTWr2+k5S125hrG+rs2nwXcXLm/U/f9X9/GnCuBxSc9KuiaXnRYRe/Ly34HT6onWMccw1OHQtC1J84FzgWcouc5mNFVl3SQ9CbyvzaqRiPh13mYE+C+wbphytbEH+EBE7Jd0HvCQpI8OQa7KTZcT+ClwC+kX3S3A90kdvR3pgojYLem9wBOSthRXRkRIqv3Z8WHJkQ1N25J0AnA/cF1EHCjeFiyjzo7pjiAiPj/deklfBb4ILM+XVQB2A2cWNjsjlzFN+UBzddhnEpjMy89K2g6cnTOcUdi00lxUUF9T9ZpT0s+AR/K30+WsQt3HP0JE7M5f90l6kHQpY6+keRGxJ18+2FdTvE45aq3DiNjbWq6zbUl6B6kTWBcRD+TiUuts1l4aknQR8B3gkoh4o7DqYeAKSXMkLQDOIt2M/QtwlqQFko4HrsjbVpX3VEnH5eWFOdeOfDp4QNIypY8FVwJVfnofqvqacv3zUqD11EennFWptf0USXqXpBNby6QHJ8ZznqvyZldRbTsq6pTjYeDK/CTMMuD1wuWQ0g1D28rv8buAzRFxW2FVuXVWxp3vYXiRbui8AmzKrzsK60ZId/63UngCh3QH/qW8bqSkXJeSruNNAnuBx3L5ZcCLOetzwJcK+ywlNcrtwI/IfxFeRa6666tNzp8DY8BofhPM65azwjZXeX10yLGQ9JTLC7lNjeTyk4HfA9uAJ4G5FWT5Bemy55u5fV3dKQfpyZcf5/obo/D0WkW5am9bwAWkS1Ojhd9dK8quMw8xYWbWcLP20pCZmfXGHYGZWcO5IzAzazh3BGZmDeeOwMys4dwRmJk1nDsCazxJJ0n6ZuH790v6VQnHaQ1zfPM02yzKQyAfGvTxzTrx3xFY4+XBvR6JiMUlH2c1cCgivtfDtoci4oQy85i1+IzADG4FWp/E1ypNEjQOabwqSQ/lyUB2SrpW0vV5lNgNkubm7RZJ+l0e7fNPkj7c7aCSPqv/T4LyfGtYCLOqHdODzpkNyCpgcUR8At46QyhaTBoO+J2koUtujIhzJf2ANPbTD4E7gW9ExDZJnwR+Anyuy3FvAL4VEevzaJP/HtC/x6wv7gjMunsq0iQhByW9Dvwml48BH8u/xD8F/LIwXPCcHn7ueuA2SeuAByJiV7cdzMrgjsCsu8nC8uHC94dJ76G3Af9snVH0KiJulfRb0qBi6yVdGBFbuu1nNmi+R2AGB0nzwx6ViDgAvCzpcnhrQvGPd9tP0qKIGIuINaRhrLveVzArgzsCa7yI2E/6RD4uae1R/pivAFdLag393Mt8xdflY46ShkN+tNsOZmXw46NmFfHjozasfEZgVp1DwDW9/EEZaXIgs0r4jMDMrOF8RmBm1nDuCMzMGs4dgZlZw7kjMDNruP8BfSS/xz+vrAQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# find the segments of good data for both source and receiver\n",
    "bb=np.intersect1d(sou_ind,rec_ind)\n",
    "if len(bb)==0:raise ValueError('Abort! no good data in overlap')\n",
    "\n",
    "# do cross correlation\n",
    "corr_day,t_corr,n_corr = noise_module.correlate(source_white[bb,:Nfft2],receiver_white[bb,:Nfft2],fc_para,Nfft,dataS_t)\n",
    "\n",
    "# plot the waveform\n",
    "tvec = np.arange(-maxlag,maxlag+dt,dt)\n",
    "plt.figure()\n",
    "plt.plot(tvec,corr_day)\n",
    "plt.xlabel('time [s]')\n",
    "plt.title('cross correlation function between %s and %s'%(ssta,rsta))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Save cross correlation data into ASDF file\n",
    "\n",
    "Though we only have one station pair, we can try to save it into the ASDF file. We save the cross correlation data into the auxiliary structure of ASDF, which has two dimentions (data_type and path). In this example, we use the station and network name of the source and receiver station to define the $data\\_type$ and use the channel names to define the $path$. The two tags are chose because the two-dimention variable are enough to define any cross component of the cross correlation functions for any station pairs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "cc_h5 = 'cc_example_TA.h5'         \n",
    "with pyasdf.ASDFDataSet(cc_h5,mpi=False,mode='w') as ccf_ds:\n",
    "    # location info \n",
    "    coor = {'lonS':slon,'latS':slat,'lonR':rlon,'latR':rlat}\n",
    "    # cross component\n",
    "    comp = tr_source[0].stats.channel[-1]+tr_receiver[0].stats.channel[-1]\n",
    "    # parameters to be saved into ASDF\n",
    "    parameters = noise_module.cc_parameters(fc_para,coor,t_corr,n_corr,comp)\n",
    "\n",
    "    # data_type name as source-receiver pair\n",
    "    data_type = tr_source[0].stats.network+'.'+tr_source[0].stats.station+'_'+tr_receiver[0].stats.network+'.'+tr_receiver[0].stats.station\n",
    "    # path name as cross component\n",
    "    path = comp\n",
    "    # command to save data and parameters into asdf structure\n",
    "    ccf_ds.add_auxiliary_data(data=corr_day, data_type=data_type, path=path, parameters=parameters)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Read the ASDF file\n",
    "\n",
    "Finally, we want to read the cross correlation function we just saved. To retrive the data, we simply need the two tags we just created for the auxiliary structure in ASDF, which are $data\\_type$ and $path$. Note that we do not necessarily need to know the two parameters beforehand, because we can simply get the two parameters from reading the file. You will see how we do it from the codes below.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['TA.K62A_TA.K63A'] ['ZZ']\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD5CAYAAAAqaDI/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deZhcVZn/P2+tvae7k05n6YROQgIkrEmzyiICAXQkbCrqYBxhcFR+o4MbDI6gOIC7ojKYARzEYQBBMQoSkrDv2UhCgOwJWTpJL+l9ra7z++Mudav6VneH6iVJvZ/n6aerbp2qOn373vM973LeI8YYFEVRlOwlMNIdUBRFUUYWFQJFUZQsR4VAURQly1EhUBRFyXJUCBRFUbIcFQJFUZQsJzQYHyIiFwG/BILAvcaYO1NejwK/B+YAdcCnjDHb7NeOB34LFAFx4GRjTEdf3zdmzBhTWVk5GF1XFEXJGlasWFFrjClLPZ6xEIhIEPgNcAGwE1gmIguNMe94ml0D7DfGHCkiVwE/BD4lIiHgD8DVxpjVIjIa6O7vOysrK1m+fHmmXVcURckqRGS73/HBcA2dAmwyxmwxxnQBDwPzUtrMAx6wHz8GnCciAswF1hhjVgMYY+qMMT2D0CdFURRlgAyGEEwEdnie77SP+bYxxsSARmA0MAMwIrJIRFaKyLcGoT+KoijKATAoMYIMv/9M4GSgDVgqIiuMMUtTG4rIdcB1AJMnTx7WTiqKohzODIZFsAuY5HleYR/zbWPHBUZhBY13Ai8aY2qNMW3AU8Bsvy8xxiwwxlQZY6rKynrFOhRFUZQPyGAIwTJguohMEZEIcBWwMKXNQmC+/fhK4FljVbtbBBwnInm2QJwDvIOiKIoybGTsGjLGxETkeqxBPQjcb4xZJyLfB5YbYxYC9wEPisgmoB5LLDDG7BeRn2GJiQGeMsY8mWmfFEVRlIEjh2IZ6qqqKqPpo4qiKAeGHYOtSj2uK4sVpR921Lfx9Nt7RrobijJkqBAoSj98/dHV/MsfVrC1tnWku6IoQ4IKgaIMkN+/tm2ku6AoQ4IKgaIMkE37Wka6C4oyJKgQKEo/1LR0AtATP/QSKxRlIKgQKEo/1DSrECiHNyoEitIHbV0xWjpjAMQPwVRrRRkIKgSK0ge1zV3u45haBMphigqBovRBTYu1R1JAIK5CoBymjHT1UUU5KGls76ahrYv6VmufpDEFUXrUNaQcpqhFoCg+XPLrlznnx8/T3RMHIC8SxH6oKIcdKgSK4sP2ujYAVwiioSA9cVUC5fBEhUBR+qC7x3IH5YQDmj6qHLaoEChKH8QciyAcRHVAOVxRIVCUPuhyXUNqESiHLyoEitIH7V09AOSEgyoEymHLoAiBiFwkIutFZJOI3OjzelREHrFff0NEKlNenywiLSLyjcHoj6IMFq22EKhFoBzOZCwEIhIEfgNcDMwEPi0iM1OaXQPsN8YcCfwc+GHK6z8D/p5pXxRlsGnvsspL5ISDuo5AOWwZDIvgFGCTMWaLMaYLeBiYl9JmHvCA/fgx4DwREQARuRTYCqwbhL4oyqDSphaBkgUMhhBMBHZ4nu+0j/m2McbEgEZgtIgUAN8GvjcI/VCUQae9q4dQQAgFRIVAOWwZ6WDxrcDPjTH97vghIteJyHIRWV5TUzP0PVMULIsgFBSCgYDWGlIOWwaj1tAuYJLneYV9zK/NThEJAaOAOuBU4EoR+RFQDMRFpMMY8+vULzHGLAAWAFRVVekdqQwLbd09hAMBggE0RqActgyGRbAMmC4iU0QkAlwFLExpsxCYbz++EnjWWJxljKk0xlQCvwBu9xMBRRlsttW2Unnjk6ze0eD7eigggBUsDocCBAIy4DLU22pb+drDq+iKaUkK5dAgYyGwff7XA4uAd4FHjTHrROT7InKJ3ew+rJjAJuAGoFeKqaIMJwtX7wbgqberfV8P2kLQZscIgiIDdg196/E1PPHWbpZvrx+czirKEDMoZaiNMU8BT6Uc+67ncQfwiX4+49bB6IuiDITqRmufgfLCHN/Xw8EAnbE47V09hIMBK1g8QNeQYwlEQyMdglOUgaFXqpKVVDe2A+n3IfZaBOGgEAgIxgxscxpHCMJBvb2UQwO9UpWspLrBsgic/YhTCblCECMUDBC0lr0MyCpw6hNpuqlyqKA7lClZyY791n4DremEIGgHi7st11DAFoaeuCEc9P/Mny3ewLpdje4eBk4Ja0U52FEhULKOeNy4K4Zbu9JZBJax3N1jCAfFtRDifVgEdy3dCMD4UVbcIaZbmimHCOoaUrKOTk9aZ0tnj28bJ0YAlq8/6LEI+sOxCLpUCJRDBBUCJeto81gBaV1DHiEIBYSADFwIHKFR15ByqKBCoGQd7d0JKyBdsHgwLIJutQiUQwQVAiXr6PAKQcdAhEASQjCQrKGYCoFyaKFCoGQd7V3WAJ0fCaYNFnvXAIQO0CJwmmiJCeVQQYVAyTqcGEFZYTRtjMBrEUS86wgOYG2AxgiUQwUVAiXrcGIEYwqi/S4oA2tNgbOOIH4Ak3x1DSmHCioEStbhxAjKCqN0dMd98/2dBWVgrSkIHUCMwEGFQDlUUCFQsg5nMdmYgiiQ2KDei+BxDYXEs7J44IO7riNQDhVUCJSso91jEYB/CqkhMfMPBbwxgoF/T3dsYNaD0Q1vlBFGhUDJOtpTLQIfIfDGhK11BNbjAwsW968an7znNU67Y+mAP1NRhgKtNaRkHR1usDgCpFlUliQE1p7F0HetoVQGIgRvbrM2rzHGICL9tFaUoUEtAiXrcHYdK86zhMDfIvC4hoLiWgQD3a4SDixGsNveKEdRRoJBEQIRuUhE1ovIJhHptQ2liERF5BH79TdEpNI+foGIrBCRtfbvjwxGfxSlL9q7e8iNBCmIWgaxnxB4h/twMHBAtYYcYgewjmD9nqYBt1WUwSZjIRCRIPAb4GJgJvBpEZmZ0uwaYL8x5kjg58AP7eO1wMeNMcdhbW7/YKb9UZT+6OjuITecEAK/CqRei8DaqnJoXEMRezvL9/Y0D/hzFWWwGQyL4BRgkzFmizGmC3gYmJfSZh7wgP34MeA8ERFjzCpjzG77+DogV0Sig9AnRUlLW5dlEeRHrR1mfC0Cz3gfCgiBDxAsHohryMlGqm3uGvDnKspgMxhCMBHY4Xm+0z7m28YYEwMagdEpba4AVhpjOv2+RESuE5HlIrK8pqZmELqtZCvtXZZFkO9aBH27hirH5A9JiYmeuHFTWXXxmTKSHBTBYhGZheUu+mK6NsaYBcaYKmNMVVlZ2fB1TjnscGIE0ZC1YthXCIzhw0eVse3Oj3HhrHEHVHTOobufonPegncqBMpIMhhCsAuY5HleYR/zbSMiIWAUUGc/rwD+DHzOGLN5EPqjKH3iWAQiQn40lNY1FPCkcx5IGWqA3HCw38G9zRObSFepdEd9G197eBWN7d0D+l5F+SAMhhAsA6aLyBQRiQBXAQtT2izECgYDXAk8a4wxIlIMPAncaIx5ZRD6oij90m4HiwEKoiFfiyBuDN6sflcIBpgJlB8N9hsj8H5vurZPrq3mibd2c/uT7w7oexXlg5CxENg+/+uBRcC7wKPGmHUi8n0RucRudh8wWkQ2ATcATorp9cCRwHdF5C37Z2ymfVKUvnBcQ2AJQTqLwLvAy00fHaBFkBcJ9WsReL83XduIvYDhybXVA/peRfkgDMrKYmPMU8BTKce+63ncAXzC530/AH4wGH1QlIHS0ZWwCPKjQVrTpI96F/oG3TLU/kKQWi8oLxLsN1icHCPwb+tYDRpDUIaSgyJYrCipPLFqF5//3ZtD8tltHosgPxqiOc2eBJ4tCfotQ50aRLaEoD+LwBKgSDCQNkbQ3NHt+/mKMphorSHloORrj7wFDE0NnnavRRAJscenvIMVI/C4hvrJGkotPZEfDVHT7JsJ7eK4hkblhdPGCByLIBY3xHrihII6d1MGH72qlIOawa7pH48bOmNx1yIIBcW3fpAVI0g8728dwQeyCGzXUGleJG3bpo6EtdLhsRq6YnEeeHWbuoyUQUGFQDmoaffZNMahurGdu5/fdED1/J0FXI5FEA4GiPlsNhM3xj99dIAWQU54ADECe7ZfnBdOO6C3eIWgO3EuHnpjO7csXMfvX9ve53coykBQIVAOatr6EIKvPvwWP3p6Pat3NjLv1y/z9q7Gfj/PFQLHIgiIb3E4A+AXLB5gjKAvv7/bly7r9aLccNq23hRTrxA4u6rt3N/W53coykBQIVAOatq6/AO5APtbrfo8z723j9U7G/nh0+/1+3mOheFYBKFgwH/mnmZBWboy1KlWRX401GffATpjVjnsvqwHJ1gM0NGd+I4m+3hDmy4064/2rh7/PScUFxUC5aCmL4vAGaerG9uBxGDdF6kWQTgoaV1D3k9zRCFd+miqRVCUY2UjpWsPlp8/EgoQDkp6i6AjxqjcMJBsEexusALcm2ta0n6+YnH2j5/j2FsWjXQ3DmpUCJSDmj6FwB6qt9dZ7hGnVHRf9LIIAoG0rqGAj2sobYwg5TOKcsMYAy19WAWdsTjRUIBoKJA2RtDcEXN3UuuMeYXAEr+Ne1uGfc/jvU0dNLQdOtVS+8veUlQIlIOcvoLFjr9+4z5rVhwOJkbuzliP6zpK+rxewWLxHYTjKWmriVpD/n3pbRFYs/jGtm530E4lYREEfLOj4nFDS1fM3VvZ6xpyPrO9u8f9m4aLU29fypk/fG5Yv3MkeXNrPS9uOLwrHqsQKAcd3hmu1yJ48LVtbNyb2MClwS7EVm8P+F6f/j/9bhkn3ba412f3ChYPNH3UtQj8Z+6pn1GUay3Rufv5zZxx57O8W917B7KunjjRUJBwMOBbqbStuwdjoKzQEQKr721dMfY2dTA6395zuWP4/d+His89NgjptZ/87Wt87v6hWdyYSlNHN79+duOwpwWrEGQJ6/c089sXNrOj/uDPMvEO/k7AtaO7h//4yzo+teB1wJot16fM+L0lG17dXOf72a5rKJJwDfXETS/3ijEkLShLrCPw73M6i+CxFdZWHXuaei9a64z1uBaBX7DYCRQ7FkGnLRZPv72HuIGLjxsHJK81GGoGY2AdTup8rMKDmQUvbOEnz2zg8RU7h/V7VQiGgK5YnF8s2TAks6aa5k5+tXRjn0HIR5fv6DXg/+SZ9dzx9/e454WBVfpu7YyxfoS2T/QWY3Nm8I4rxHltf1tXr8G32WdATJ1ZpcYIHHdS6kBsjEkTI0hnESQfL7IDvM7nRnxWBHfZMYJIyHINpYqRM9N3YgSORfD4yp1MLs3j3KOs+ozDOTtPFd+Rpidu+MYfV6dNHfauGu8vndePdPGXp9/eQ+WNT/q6HwFWvb9/QOnMqThbl26pbT3g92aCCsEQ8McVO/jFko38doCD7oFw05/W8NPFG1i1Yz9g+S/vfWmL+3p3T5xvPbaGy+5Orupd19Jp/x7YjXzdg8u58BcvJgUoD5T/eWUrD762zfe1upbOpCwYL96BzbEOdtlC4OwzvM8OADoZNQBNPjX7U/3nqa6hoB1gTh3I4ymuoXBQEEk/mKSzCFL/Di+ddowgkkaMnBpICddQnN0N7by6uY7LZ0+k0P6OgbiG/rxqJ29ure+3nR+b9jXzyqZaAGpaEoHXoap/tH5PM79csnFAbbfUtPDYip189eFVvq/v9VhifcWb0uG1tryicPfzm6y+7vWfLF1296v8w69ePuDvywlb1+OuNHGloUKFYAjYa89CBrdCjsWOeusC2dvUSXdPnE/+9jV+8OS77k3pbHZSmzLgO/nm9QPM9nhlk+VaqW5Idmm8saWOuT9/od8ceYBb//oO//GXdb2Od3T3MOcHS/jWY2t83+etBuoMoI5FUJBjCcH7tsVz3dlT3baOReC1llJv/nQWQaqP35Bca0hEiIYCSWUevKSLEST+jt7nqzMWJxK0XEPQ23ppdi2CRIzgr6t3YwxcMbvCFcWWzr7XErxf18a/PbI6qYjfX1fvHnA2zfk/e5HP3vsGkHxdpWYOdcXig+LbvvAXL/LzJRuS1lB4+dniDfxtzW7iccPmGmvmHE6xuD694HU+dtdL7PX8jW3dB245eYWkuTPGJjsxwbnf9vq4/Pz4zhNreWlj/wHnFvva3+mx6O9/eWvSZG8oyEoh2Li3mT+vOnAfXGN7N19/dLXrMqlr6WTh6t29zEdnN6ncSPqafnc/v4nrfr+cLTUt/HLJRr768Cr2NHb068N3LuZtda28sD5xYTk3dWuaAdoRgANN+9uRsnL1R4vWs2FvC2+93wDA9rpWXt/i74936OjuSZo9Pv32HgCeW78vqd2Ta6r5zhNrafYMbO3237PLFiRnAN9mm86fOWWy29ZZZOX1C/cSgu7U9FFbCHq5hiA1GzUnHExrxaTOjp1BOl0/wE4fDQddd0DqIOrM9F2LINbDtro2xhREmVSaR6Etiqkxgi01LfzHE2+7fXrw9W0AjC2M0hM3vF/Xxv/7v1X8yx9W+P4t6ejo7qHWM7B63UTxuGHODxYzP01Qdf79b3LrwnW8uqmWB17dlvY73tuTCKq3dMZY/M7epLhET9xw19KNXP/QKm796zr3b8izLTyH17bUsW53E42e672vVGSA376wmUeWvZ90zDvQ3/T4Ws7/2Qvsbepwr5ed+xMzd2MMD7/5fi9xaGjr4g+vv8/V9/mfm7d2NHD+z16gqaPb/Z9v2NtC5Y1P8uBr2/j+397hB0O8MVFWVh+d+4sXMQbOO6a8lwnfF9/842qeeWcvRbkhbrz4aOb8YAkAE4tzmHNEqduu1r5BHH92Y1s36/c2c8qURJt7X9pKfWsXJ0wq5r6Xt9DUEeMvb+0GYM2tc3371RWLuz7P7bVt5IUTF/+uhjbGjcrpNfPc1dDOtQ8sT1gErX3PHn+yaH1SGYUd9e00dXS7/Rk/Ksc6vr8NYwzn/+wFunsMW+/4aNoqoTf/+W0eX7mTFd85n5K8iCsER44tSGr3lYdWAnD61DHuMdc1tD+RLgmwra6N0vwIJfkR3vruBfz3S1v4zXOb6e6JJ92IqTd/e3cPkWDAreLp/E4NgsZTa0wAOaH0QpAqJKFgwE5NNb79ADt91GMRpKaQOjPi0Z700f2tXZTkWf8LRwhSXUPXP7SKd6qb+MfTjuCocYVss9dZ5ISDfPXhVfxtjbXJzYHGgHY1tFPb4i8Ef12zm+aOmG+QPh43vLa5jhc21PA/tgjMP6PS9zveq070aem7+/jOE2/zxXOmctPFxwAJSxBIqrOUrjjhLo9FW9tsuSPjcRhTGGH8qNyktnf83VqZ/qmTE5MLb4zBmbjsbmhnX7N1/MeL1jOuKIcr5lTwTnUTN/5pLRfOKk/0KxZPm7jgcPuT77JpnzW5cqw75zr/6eINbrvOWA/RUND3MzJlUCwCEblIRNaLyCYRudHn9aiIPGK//oaIVHpeu8k+vl5ELhyM/vRFbUsnzji3Yvt+3zY76tvcQby1M0ZXLE48bnjRNu027m1xFzEBVNsXS2N7N8u21btmXWN7N4vf2csJ33+GT/72NdfvZ4xxP39PY0evFbEPv/k+R/77U7y1oyHp+NbaVndg2VrX6vrJwZpB/PSZ9Ul72xpjeOiN7W7qYm44SENbV9oAWDxu+PVzm7j7+URs47EVOzj+1mdcH7Hjk//9a9uZdcsitz9e33Eqj6+0rK85P1jCHX9/172Z08Ur1u1OBNnaUmrqOHGA7XWtHDE6D4DivAhl9mDZ0hFLunnbU9wB7V09rh8WPMHiXv7u5GAxWP5bby6/Fz9/udfn75fr3xXrIRoOuIHk1PiDEyspzAkRDQXo7O5hf1sXJXbaaL7rGor5vs+5Fpzr5L09za4IQGIdxqZ9zfzr/63q5Sr625rdSf+L3Q3t7rUOVsAerJXdC+1JzMTi5MG1qaObx1bspKsnnjRrTxXU7p449a1dSfeV02+v5ZtOvPY1JfrudSmt2Zm4hz614HU+dtfLfPzXL3P6Hc/y4Ovb+eYfVwPJ/v+XN9by7cfWsG53o7tGxcvmmlb2e0p7fN3+jOXbrPHEey3c/tS7fPl/V/r2zcGJwzV1dPf6XwY9k6utQxhAztgiEJEg8BvgAmAnsExEFhpj3vE0uwbYb4w5UkSuAn4IfEpEZmLtcTwLmAAsEZEZxpghWyHjDZgt21rPuUeNZVttK198cAVzZ5Uz/4xKzvrRc8w7cQK/vOokZt2yiDOmjeaOy4+joztOKCC8vKmWy+9+1f2cbbWttHbG+OcHlvPmtnrX7dDY3s0Ndl19gK01rUwszqW1q8dNBdzV0O7mwzvc8ff3MAZ+tXQj933+5MT7a62LckZ5ATvr25hcmkdRToimjhj//ue1GEOSa2l/m3UTOkwZk8871U00d8YoygnzxpY6vvPE2zz+5TPYVtvKI8t2JPWjOC/MStsFtHZXIx86cox7oa7bnZwXv+jtPTy8bAefrJrE/DMq0wZV//ulre7jmuZOjDFsrmlJcqOttbMtCj31epyboKkjhjGG7XVtSRZWgW2xNHfEkmaqvSyCrsSmNJBYjexnEaQaOH25hvzKVHzl3GmU5EW44+/vpY0RRIMBwqE0wWJ7pp8fCbnfvb+tiylj8gHLL54bDvYaXJy/xZm11qaJBXR097C1tpWr73uT6sYOqipLmHfCRL780Apum3cs1z+UHIB96/0GHluxk+MmjmLtrkbqWrvoisU5/Y5n3TZJ+zDH4nzojmfdoPftlx3HwtW7efa9fdS2dFJRkue2veOp97j/la2ce1SZe8yJC22pbaUnbpX82JASnC3OC9Mdi1Pb0unu1+CdCKRep17+44m3AZg1oYjTpyWs0H+8z4qHPLJ8h+/73tqRPIGsHJ3HYyt2cstCKx7mXc/yzLo9Ke9t4KzpZUnHnLGguqGjV+ab1zrftK+Fo8cVpf17MmEwLIJTgE3GmC3GmC7gYWBeSpt5wAP248eA88TyI8wDHjbGdBpjtgKb7M8bMpwZTeXoPJZtq2dHfRtXLXidzTUt/OrZTXzJ9jm+taPBvahf3VzHO/YFdfnsiUDyBf+TZzZwzo+f481tlsg4s7/tda00d8a49swpgDXTPfX2JfxpZWJwXr+nmdQJuvN8474WfrzoPXdV49Zaa5Cfc0QJ+9u62dfcSeWYfIrzwu57vDOVh5e9z17PTGlqmTWANNjuoa88tIqN+1pYvq2eS379Cv/7RrJ/tMrj7nIGQOdCPWv6mKS2v3tlG+t2N3HLwnW8trnON6XxzCMT7ynNj9De3UNrVw+X/uZVPnRnYjBZs9MSgjGFUdq6emjusP7W0vwIPXFDbUsXuxra3QEREj7/jlhPUkDcL0aQ63GphfpMH01Wgmg4mDZY7GcRfPPCo7n2rKnkhYNpXUOWRRC0+9A7WFwQDREMiGuN7G/rptS2CMAKnqeea8e62ddkCW26oHDcwD/e+wYtHTFK8sK8ubWe17bU8cqmOm54dHWv9n94YzstnTH+87JjAavoX11r4rNnlBfQ0hnjsRU7mXLTk7y8qSZp97ejxxfyj6dZbpfUPv3+tW0APLe+xrU6nUlNVyzOtH9/il8u3ciWlNpKf/nKh7jpo8cQN4nY0G6fjYb64ta/vsOvnk2fpeRYOc6A7QSMv3fJLKaOyacrFufvnj2ldzcm4gZed2lxXpgHXrXcWW9ureeyu19hR32b2+9dDe20dMYoyklMirzC4CRwDAWDIQQTAa907rSP+baxN7tvBEYP8L2Dyr6mDiKhABfMLGf1jkZ+++Jm9rd18cAXLP1ZZpt34WCArTUJU+zel62Z7DcvPJqvnDvNPT7B9pmnZukAvG2Lx8lTSgkHhRc21LC3qZN77VnxuKIc113kVzDt/fo2fvPcZj53/5vE44attS2UFUaZWJxLe3cPu/a3MdZ+7uDNwb/n+c1uDjrAVHvg/Nva3XT3xF1/74a9iZvLGSTPnlHGNy88yj3uzLKa2rs5fepoHrzmVKaVJQbiLbWtREMBJpfmcf1DK5NMcoAr51Tw+y8kNP7EScXu56YOZI5LY1xRDo3t3Wyx/w/Oe1a9b/2PZpQnYgyOu6e9K7m0hF/6qNf6cPzzvumjJJMTCvRhEaRPpcyNBNMGiyN2LAH8XEPdbtA5Jxykwy6bUZKX+J8WRkO9ZpHOd+1r7qSxvZuunniSaN5+2XF86cPWNbyroZ3zZ5Zz1vQy3txa714/qW5JsDLVAgLHjC+iIBqirrXLHdA/d/oRXHrSRHrihu/9dR3GkGSNAkwozqWswLpfapo7aeuKsbW2lT2NHUnnz1kx7QRijyovBOCPy3ckuSCL88JMLs1jrB1Md9xDexqTUy9He4QzHX2l1p4waRSQEHsnU2l6eQEXzCqntqWLzTUtXHzsOI4eV+gWBATr/E4dk8/D153G58+oZMm7e1m2rZ57XtjMqvcb+PL/rnTPYXVjOy0dsaT/lXNeSvLCPLLsfTbtG5q1PYdM1pCIXCciy0VkeU3NB6/7sbepg/KiKCdXltLVE+ehN95n5oQizpg2mkJPpseO+jY21SRO+ort+zlidB5lhVFuuCAxQKa7/3PCAffCOWJ0HpNK8njDvtgcH/msCQkzb4Z9safjR4vWs6WmlSlj8hllDwSba1opK8xJEgJvXZumjhhXzKlw/64p9sD9o6fX88y6vW47x9y+b34Va26dy5++fAa/+cxJHDWukJe+dS7Txxa4s6ymjm43NTInnBy4mjWhiN/908nUtXa51oWT+lheFHW3e4TEoO7NEvEiYr2noa3brbDpvGeFLQRHjk2cM9ci6O7xdQ11xeJU3vgki9/ZS64nRpA+a6j3Fpk54SCd/WQN3TZvFi9+89yk1/IifVkEQcKhdMHimJsumxsOWtkqcZMkBKkWgdfPvKexnZtt94cTmA8GhM+cOpkvetJuK0pyOb5iFPuaO9PmxTuMK8ohHAxQkh9mf2uXO5m4fHaFm1DgWDZPrd1D5eiE+6coJ8yYQqvvtS1d3PSntZz7k+f55VJrNn7f/CrKCqN89rQjAGsQLcwJ8dRXz+KaM6ewr7mTHfXtFNvB8uMrihERN2bS0G7933fubycguAJdHAwAACAASURBVDEkJ+sqlc+emggK7+sjlTY17uEM3HmREGUFUbp64myra+PIsQUUREO9Ft1dPnsip00dzT+fNZWJxbl844+rec0OIK/1LDqrtidFR4zOJ5X/+IeZxA0s8ty3g8lgCMEuYJLneYV9zLeNiISAUUDdAN8LgDFmgTGmyhhTVVZW5tekT4wxPL5iJ0+8tZvywhxOrrTcHnFjDWAi4t6QJ0wqpjMW5/XN9QQEfvKJEyjMCfGjK44HrJvpjsuP47dXz3EvvlSmlSVmqxUleUwqzevlPvAKwdHjrEHthApr9uEdgy6fPZF7XtjM8u37mT62IGkR1djCKBNLPEJgD9jzTpwAwCerJvG5M6wba/bkErfdox7/pyMExXlhwsEAsyeXuIuVJpXmMWVMPi9vrOHtXY00d8TcG36S7eP9sO3XLcmLMK2sgNH5ETbb5rOTZZRjZzucPnU0ACdNtgb1tTuTV186f3duOEhxXoSGti7e2tFAJBjguInWuVmxbT/hoLg3OlhuG7Bm/PvbulxLzZkdL9+WmPF5YwTpcvhTaw0B1jqCNMFiZ+Z2+rQxTPb0y/q+UJoFZVYGkxMsTq031NIZczODyotyeNfOqCnxzHBH5YaT9iTwBlOXvLuPJ+3g8HRbCBzB9F5DFSW5TLAHOz9LwNveudZK86NJFsGYgojbV+85Or6imOvOnso5M6xrZHS+NSi/X9/mZsn935vWpOGkySUsu/l85p9+hPv+krwIwYBw6pRSYnHD+/VtzBxv3TfOvVJs9805D+9WNzO1rMBNIHAmI14evOYUvnzukb5/K8AZ00Zz7ZlT+OLZU/lKmnb5kWCSyEwrK3CF28tke2DPj4b49sVHs72ujfbuHj5y9Fi3zdHjCtla20pTezejCyJJk1KA6WMLmTWhiBeGqPjdYKSPLgOmi8gUrEH8KuAzKW0WAvOB14ArgWeNMUZEFgIPicjPsILF04Ehqe4kItz/iuWSKR+VQ0l+hKPHFfLenmZ30HaCbHNnlrN6RwPPrd9HRUkeV86p4LKTJia5bz5t56+nGxiOHFvAut1NlOSFKYiGmD62gBc21CBiZQLE4oZTpozGCoskLILjK4pZvbMRY+C2S4+lqb2ba86cwp9WWvo4/4zKpGDY2KKoewN6ufljx3D57AqmlRXw9QuO4ovnTKMoJ8zb37uQ425dxAsbLF9sbjjouoZG5fqb0HmRIHED//Crl8mPBF2RuOPy45hzRAknTi7m+fU17gBbXpTDO3am0uTSPNbuaiRouz/u+cc5vL61juMrHDePNfDcftlxPL5yJ7nhIC9vqrWFIExTR4w/r9zFBbPKGW27udbsaqRydH7SIiLHNeT40SeW5LK7scN1DXnXLOSGE+crlHZBWXKtIes7LPeMH07pCT8XX14k2Ct7KdYTJ25wS0yAf7DY+d9OKM51BwEnfRSsc71xb637/IlVu8gJB5hzREmST/koe6Lh/I+81s6kkjxy7ONv7WhISnt1OMn+Hzt9Lc0L89z6Gl7aaH33mIJor7UTYM26T7XFH6wSCsV5YR628/UnFueyq6Gd3HDQ/btCwYAbE3Fm/yfY1iBA1RElHDm2gCtmVwAJkWps7+bJNdU8t34fHztuvBug97pHRSyRL86NuC4lP37+qRMpL8pxn/tZdbmRoCs2YAuBzznwujA/dtx4nlpTTVVlCSV5EZ59z7ouP37CBH68aD1guftGF0SS4ivFeWHOnlHGf7+4heaObvceHCwytghsn//1wCLgXeBRY8w6Efm+iFxiN7sPGC0im4AbgBvt964DHgXeAZ4GvjKUGUMXzLTye52UrH//qJWb7AQxv3iO5Teda7fb19zpztTTbXryh2tO5WPHje913DEnj7MHvBPtGbAxcHJlKfmRIGdMS9wgzuzm1KmJAO3Vpx3BV849kpxwkHs/V8X3LpnFjPLCpNncxOJcKkqSTddIKMDYwhx3FhYIiDuLL4iGGGdf4MdOLKK8KOr6pr2f62XWhFHu49auHtc1VJIf4Z/PnkrVESXc+vGZfPfjMwEYNypxA33l3CP54jlTmX96pfUdeWEunDWOUblhCnNCrLTdPOcdM5bHv3QGY4usGys3EnRnes2dMa70uB66YnFXFBycma5TfnpsUQ7BgNDWFaOxvdsVUuezHZz/a2+LIF36aN/rCEJphCB1EHHcQE7ROb8+WDe8da4nFifOabHHNTShOJe9zR1091irev+2ppoLZ41LykxZcPUc1/rMDffOQ68oyXMtt65Y3LWWre+1ri1HuJ2JT2l+8iCaEw4mDU6XnjiBH11xfJIIOFSOzqehrZsxBVGunFNhf14kSZwKouGkv9U7aI8blcv35x1Lpe1Ld+o6balp5SsPraQnbjhmfJF73Nsv5xoalRvutRrZsTiBpAHeex685EVCjPH0a+aEIt8B+iiP2zcYEO65eg7XnjU1SWicJBSw3H2fqJqUlBVXnBfmitkT+fVnTurV78FgUBaUGWOeAp5KOfZdz+MO4BNp3vufwH8ORj/649yjxvKLJRvdE3n2jDI23/5RdzD48oencd3ZUwkFxL15jxnfd7rWmdPHcOb0Mcx8bhPBgHCnvSjl0pMm0twR45sXWfGEEz0zmq/PncH6vc0EAsK/fuRInnlnL6dPG83f/t+ZrmVQmeJeOH9mYpFKcV6yWe/sfeuQH+l70cmUMflUN3bwqZMn89fVu7Fi9+mF4J8+VEl+NMS//3kt0LuOjojw+Q9NcZ97L/CxRVF3MVAqFSV5vFvdREAS5rtjEjuuIYdjxhe5s1G/vjrxivYuK2todH6EXDtb5xdLNlDf1mXl4sfiSTEC51pIddulTx/tO2vIb8KQGw72ypLptD8nGkq4hlJFpq2rh7xIwiJwmO6ZYU4YlYMxVuzr/fo2Gtu7ufjY8cyaUORei3NnjaOju4dJpbl8b94s973HjC/i3eomxo3KISDWRjxxQ9JCK2eG67hjnPTO0vze14p3NnztWVM51jOwevnyh6dx3YMrOPeoMjeTLfX8F+aEqG3pdK0EESE3HKS9u6eXzz8nHCQnHHDX+YCVoeT68qPWtfHZUyfz8qZaGtu73evnxouPJi8S5P26Ns6aUeaujA6k/B+PGJ3Hxn0t7jkCS+Anl+Zx6pRS/vW86QQD4gp3xN5jIhSQtAstx41K/B3jR+Vy1vQxbNrXwqlTRnPCpGLOmDaay+5+lVBAKIiGKBxbmBQXG0yyamXx8RWjuOvTJyWlMXpvXBFxMzgGKgQOjh/RuflmlBdy26XHuq87M4ppZflUVZZSZc+6bph7FDfMtcTCuXHum1/FzAnpv9c7CE4oznUHFYe8PkpbgBXz2NPUwezJJa7vPNdT6iCVUDCQ5M/0c0V5GecRAj9T2aGiJJd3q5soK4y6/wdnkVRuJMgoj+CNLYzS48mzTRUjZ6bb0mlZAMV5EXIjQdbvaeaNrfV8+pTJ7G5ot1xY3vTRdMFieqeP9rWOwHFB9WcRPLW2mj8u38Edl1vxpkgoSL49UPmtgnYWYXmFwPu3j7ePVzd2sPTdfURCAc6aPob8aIjrzp7KCfZMPicc5KVvfSTp8/9wzSlsrml1/+/OADe9vIBvzJ1BdWMHZ00v41/+sILTp43mtZs+Qnmh9b/1irSD97pIXTXu5YKZ5fzXZ2dzxpFj3PUhqVlbzvVQ7LnWS/Mj7GpoT3L1OBTnRty0zvvmV/HhGWWssheMRkNBtt35MQAu+fXLiCT6+i/nJDIA0yUuQCLgPH5UrpvpFw0FEBEe+eLpbjvnej9mfCGXnTSRs2ekj2eO9dwnAA9ec2rS85Mml/DSt84lNxJMKyaDRVYJgYhwyQkTBtS26ohSnl63h2PGH5gC33bpsW59nNTvfulb5/Y7iIJV+qIvvOZnXiREXgTu/3wVL22s5XevbOu3YuiE4kRw0LFA+tvlqrwoypQx+Wyra+3zJgfczBDonVnkxbm5PuRZzOME2wIiSYIXCAgBEpZaOotgX7O1crwkL0xeJMiq9xvoiRs+ffJkFtiFu/zSR1PdMnErSJD8HbZFkZpR1NzRzW+e28SUMflJOf4O3mCxs8rUGXSioYDrqmqz/wd7Gq3Mto7uHvfvclw3qX+34zLa3dDOO7ubOHZCkSumjuszHaMLom75Ci+Xz57I2MLEIOUMot7vdjKTTq4s4dsXHQ0kC0Ff/3cR4WLbnepYvl6LGRK5+hcdm3C7zjtxAnc/v9k3C2hUbpg9TR1MH1vg3j8h19qLJ7Urygn3mvFb/U/vd3f+r6Fg8sQxFUcIxo3KSbKS/UgNCPsxqTSv3zaDQVYJwYHw408czxVzKnxTufri6tOOSPvaYP1T/dwPHzm6nHFFufzulW2+axrSMVCLR0RYcsM59MRNWsvBwYmr/MBjEfnhzJ4vOnace8y5keLGJM0GHYpywrR19bj+X4eo3ad9dp2h4ryw60oAa5ByXGZ+C8p6rQMw+C4oAyv/3zvQba9ro7ali+9dcqw7+HjJiwRp74olCbRTfyYSCrgWXHtXjHW7G/nYXS/zn5cdS0d34nsqSvL4VNUkrj49+fpy3Di7GzrY19zhBoU/CLdfdhz727qSRCAdTgD0a+fPcK3bvqy/dBTnRXj0i6f3mnDd9emTqGnu5HRPHO0bc4/ik1WTklYkO4zypJQ6uP9bj7U3riiHvUX+C86chVxH+5xDJ9vJb39rL45lM66o/3MoIlx20kRmTy7ut+1Qo0KQhsKcsBtcPlRwbqaBzDQc/C76dAQDkjZo7mXOEaWs/u7cJNeOH189bzpHji3gfI8FlG8Pij1x4+t+cGZ+qTPjQMAqE+3kgxdGw0kDU1406M7+cyPedQTpLALTe0GZZ62CVwicFaeOiyeVvEiQtu4ed8U04NZuioYCrjC1dvbwxhbLVbdye4P9nVb/ggHhh1ce3+uz86MhRuWGrUJoTZ29yhccCJ/x5NX3x6UnTuSEimKmetKkQ8EAF84q56M+yRN94Q2KOvhZ7oGAuAHiVJzrwVn8BRAO9F6fcePFRyeVOfdSmBPmV58+KSlhw+GTJ09ixfb9zJ1Vzlcffsvn3RZOGemBTvp+/qkTB9RuqFEhOET5+1fP6jUrFhFev+m8XtkufeG4EfpKpfsg9CcCYLkmPmdnEzk4rqGeuLEyi6IhbvK4OJyMpdR6/2AN1I4QFOWGk3K6C6IhVwi8M/2wz6wRrPTR3jEC6/2dKfn+zkw/XWXI3EgQYxKL/YpyQrxn5/tHQgG3hER7d4+7SthJCPDL8kll/KgcNte00NwZc7OuhhoRSRIBh99eXTUs35+KYz16LYJZEy1r1+t2stxh6T/n42lcxwXREL/57OykInx+fOFDlTS1d/PZU9N7Bg5GVAgOUdK5dLypmwPl9ZvOS6rIOZI4s3hjrFnw2u8lF6T1pv+lkhMOuK6hwpyQ+1ki1oI2t5SDZ4YY8pSY8Jb5jRvTO2solLAIvDjCEE1zDp1y4U611bLCKE12mQLHssiLhOxN6S0hcwrJDUQIJhTn8oa9J0T5ANw6hyNji6LkhANJLqYzpo3hpW+dO6h+9v7K1hfnRbj1kll9tjkYOTjufmVEGTcqx9cNMxI4FkpPaiU+G8cK8rshc8NBd5OWInudAlgDcSAgvj7jsG0+/WnlLo76ztPuhjfGr9aQ6xpKsQg8qaB+ODEAp/SAd6VroaeExOMrdrH4HauEgCMafQVdHSYU59BqB6OHyyI42Pjns6by+JfO6GWVDXawdSDJHociKgTKQYUT0I2nKeLkBPT8LYLEIOC1CBxx8Ss57VgETh2oP63a5dam711ryD/f37Ew+nINQaI6pjfrpdBeOJUfDSZlbjk1fAYiBN68//IBBCkPR4rzIkkLH4eK/A8QED8UUCFQDirclL/+LAIfIXCyekSgIBJyV6c6Lpt/ON4KYl7sCWZ60wEBlryz1y3p7begDHqn2jqF6NJbBNb76ludujw+FoFtNZwzo4wpY/LdzK+BuOy8q16z1TU0XAzFqt6DgcNT3pRDlvGjciiMhtzc9FSmlxcypiCSVIHTwVkxXBANEQiIGyw2JvFeJyfeIXUB2Pa6VhwJSg0WOzP71JLS/cUIcl0h6EIkuSyyM8N04gij8yPUtXa6geWBxAicQnAnV5YMKEivZM6lJw5sPdKhggqBclCREw72ChB7+fjx4/n48eN9F/M4M3YnflAYTRYCP0IpO9S32hvhQO8YgZPamroC2BWCNK4hJ0ZQ19pFbjiY5F5w1mQ4qael+RHyIyH3M3P7KRcCMGdyCb/69EmHXLrzoUrqZOJwQIVAOaToa6m9k9XjuFsSFkF6JQgHe3/eHjvzKHX1qePiSd12MpE+2p9ryBICv1LFjuyU5EeS1j8MJEYQCEjatEdFGQiHp8NLyUqc2bMTP0isUk7/HpHEIjmndMaeNFsduqUgUi0CO2soksZ/7Lh3Gtq6yQkHfVfgOmIyOj+SZDEMxDWkKJmiQqAcNjjlFZwUUNcioO+yAI7F4JRq3utYBCnWR1+uoUgw4Fu/BhIWAVjBX7/MEycTqTRFCNLFHRRlMNGrTDlsuObMKZx/TLlbUbJwABaB93VHCKpti6B31lAAEXoVFbQWoqW/lbzVYHMjQQp8SlEk6vxHkl5Xi0AZDlQIlMOGcDDAvfOr+MKZVtXH/AEEi7185Ohye29gK80zdX4vIuSFg+7iLQdr7+H0t5IjIGAN7E5aq5dEcbxwkkUwkBiBomRKRkIgIqUislhENtq/S9K0m2+32Sgi8+1jeSLypIi8JyLrROTOTPqiKKk4rpxZfeztAHDF7Aq+87FjOHJsAaX5EersxVypriHw33+4MxZPmzEEiU1VwBrY/YrTXTTLqsBaXpS85ePhmreuHFxkmjV0I7DUGHOniNxoP/+2t4GIlAK3AFVYtbxW2HsVdwI/McY8JyIRYKmIXGyM+XuGfVIUwCp898h1p/W5yQ/ATz95gvs4FBR3tyy/BCWnpLQXSwj6HrCdfRRy0wSL/+2CGVx71hSK8yIHTbkPJXvIdLoxD3jAfvwAcKlPmwuBxcaYemPMfmAxcJExps0Y8xyAMaYLWAlUZNgfRUni1KmjD2ijbyGxqtkvVTUv0ts11Nnd0+8eDU7GUW4ksY7AuyI4GBBXAJztIBVluMjUIig3xlTbj/cAfitaJgI7PM932sdcRKQY+Djwywz7oygZERCPReDzumUR9LBzfxtjCqLkhIMDsggcN1Vu2NoX4ZdXncicI3w9qYwuiHL0uELer2/L6G9RlIHSrxCIyBJgnM9LN3ufGGOMiAwwLJf0+SHg/4C7jDFb+mh3HXAdwOTJA99AQ1EOBBGrBDXgu69DfjRES2eMj931Mp8/o5J/u2BGUvnqdEwqzeO9Pc1u8HfeiRP7bL/w+jN7beiuKENFv64hY8z5xphjfX7+AuwVkfEA9u99Ph+xC5jkeV5hH3NYAGw0xvyin34sMMZUGWOqysrUdFaGhoCIW6bazzWUGw5S39pFY3s3a3dZm5R09pM1BHCSvR1hf/tJO0Q8exkrylCTaYxgITDffjwf+ItPm0XAXBEpsbOK5trHEJEfAKOAr2XYD0UZFAIifVoEeZGgu/J4w15rl7HO7v5dQyfYO2ett3cmU5SDiUyF4E7gAhHZCJxvP0dEqkTkXgBjTD1wG7DM/vm+MaZeRCqw3EszgZUi8paIXJthfxQlI0Q8G9n7BYujiYJwO/e3s7fJ2jS+P9fQ8RVWrfwD3c9XUYaDjILFxpg64Dyf48uBaz3P7wfuT2mzE/94nKKMGCLirjT2DRanLPA69falAP1mDRXmhNn4nxfrugDloESvSkXxEJDE7mh+C8ry0uxQtWlfS7+frSKgHKzolakoHgIirmvIb0HZmAL/xV7BNAXnFOVQQIVAUTwkWwS9Xx/r2W+4xN4N7BNzKrj7s7OHpX+KMhSoECiKBxEhFreCweITJSjz7Ak8e7K1IOwzp05mgmeVsKIcaqgQKIoHa0FZ4nEqXovg5CmlhIPC1DEFw9Q7RRkadKtKRfGQVGLCRwnKPELw+TMqOWdGmW4YrxzyqEWgKB4CgisEfjEC7/4AOeEgx4zvu7KpohwKqBAoigeRvstQK8rhiLqGFMVDQDxlqNOsdywrjLqZRYpyOKBCoCgeAiLuIJ/OInjl2x8Zxh4pytCjQqAoHry1hvyCxdB/OQlFOdTQK1pRPHgtAl0srGQLKgSK4kFE+o0RKMrhhgqBonjoL31UUQ5HVAgUxUNA00eVLESFQFE8JKWPqhIoWYIKgaIkIZg+NqZRlMORjIRAREpFZLGIbLR/l6RpN99us1FE5vu8vlBE3s6kL4oyGHjjAmoRKNlCphbBjcBSY8x0YKn9PAkRKQVuAU4FTgFu8QqGiFwO9L+9k6IMA95dyTRYrGQLmQrBPOAB+/EDwKU+bS4EFhtj6o0x+4HFwEUAIlIA3AD8IMN+KMqgEPDcEWoQKNlCpkJQboypth/vAcp92kwEdnie77SPAdwG/BRo6++LROQ6EVkuIstramoy6LKipMfrDlLXkJIt9FtiQkSWAON8XrrZ+8QYY0RkwJW4ROREYJox5t9EpLK/9saYBcACgKqqKq34pQwJXteQyoCSLfQrBMaY89O9JiJ7RWS8MaZaRMYD+3ya7QI+7HleATwPnA5Uicg2ux9jReR5Y8yHUZQRwjv4B9QiULKETF1DCwEnC2g+8BefNouAuSJSYgeJ5wKLjDH/ZYyZYIypBM4ENqgIKCNNctbQyPVDUYaTTIXgTuACEdkInG8/R0SqROReAGNMPVYsYJn98337mKIcdCRnDakSKNlBRmWojTF1wHk+x5cD13qe3w/c38fnbAOOzaQvijIYiMYIlCxEVxYrioektQOqBEqWoEKgKB7UNaRkIyoEiuLBO/arDCjZggqBonjwxggCWmNCyRJUCBTFQ0AtAiULUSFQFA8BLTGhZCEqBIriQReUKdmICoGieBDNGlKyEBUCRfGgWUNKNqJCoCgedB2Bko2oECiKB40RKNmICoGieFArQMlGVAgUxYMGi5VsRIVAUTyoa0jJRlQIFMWDd/BXi0DJFlQIFMVD8sriEeyIogwjGQmBiJSKyGIR2Wj/LknTbr7dZqOIzPccj4jIAhHZICLvicgVmfRHUTIlOUYwgh1RlGEkU4vgRmCpMWY6sNR+noSIlAK3AKcCpwC3eATjZmCfMWYGMBN4IcP+KEpGJA/+qgRKdpCpEMwDHrAfPwBc6tPmQmCxMabeGLMfWAxcZL/2BeAOAGNM3BhTm2F/FCUjAmoRKFlIpkJQboypth/vAcp92kwEdnie7wQmikix/fw2EVkpIn8UEb/3K8qw4R37QwENoSnZQb9XuogsEZG3fX7medsZYwxgDuC7Q0AF8KoxZjbwGvCTPvpxnYgsF5HlNTU1B/A1ijJwvJvRqA4o2UKovwbGmPPTvSYie0VkvDGmWkTGA/t8mu0CPux5XgE8D9QBbcCf7ON/BK7pox8LgAUAVVVVByI4ijJgvJlCahEo2UKmV/pCwMkCmg/8xafNImCuiJTYQeK5wCLbgvgrCZE4D3gnw/4oSkZ4YwRBDRIoWUKmQnAncIGIbATOt58jIlUici+AMaYeuA1YZv983z4G8G3gVhFZA1wNfD3D/ihKRnjHfhUCJVvo1zXUF8aYOqyZfOrx5cC1nuf3A/f7tNsOnJ1JHxRlMFGLQMlG1AmqKGlQIVCyBRUCRfHgtQhCKgRKlqBCoCgeNEagZCMqBIriwbuOIKhV55QsQYVAUTwkFZ1Ti0DJElQIFMWDjv1KNqJCoCgeRCuOKlmICoGieFCLQMlGVAgUxYNuT6lkIyoEiuJBdUDJRlQIFMWDWgRKNqJCoCgetPK0ko3oZa8oHjRrSMlGVAgUxYN6hpRsRIVAUTxojEDJRlQIFMWDCoGSjWQkBCJSKiKLRWSj/bskTbv5dpuNIjLfc/zTIrJWRNaIyNMiMiaT/ihKpuiCMiUbydQiuBFYaoyZDiy1nychIqXALcCpwCnALfb+xSHgl8C5xpjjgTXA9Rn2R1Eywik6p4aBkk1kKgTzgAfsxw8Al/q0uRBYbIypN8bsBxYDFwFi/+SLdfcVAbsz7I+iZIQjAOoiUrKJjPYsBsqNMdX24z1AuU+bicAOz/OdwERjTLeIfAlYC7QCG4GvZNgfRckIRwDURaRkE/1aBCKyRETe9vmZ521njDGAGegXi0gY+BJwEjAByzV0Ux/trxOR5SKyvKamZqBfoygHhCMAohaBkkX0axEYY85P95qI7BWR8caYahEZD+zzabYL+LDneQXwPHCi/fmb7c96FJ8Yg6cfC4AFAFVVVQMWHEU5EByLQGVAySYyjREsBJwsoPnAX3zaLALm2gHiEmCufWwXMFNEyux2FwDvZtgfRckIjREo2UimMYI7gUdF5BpgO/BJABGpAv7FGHOtMaZeRG4Dltnv+b4xpt5u9z3gRRHptt//+Qz7oygZoTECJRvJSAiMMXXAeT7HlwPXep7fD9zv0+4e4J5M+qAog4laBEo2oiuLFcVDQNcRKFmICoGieBDNGlKyEBUCRfGgMQIlG1EhUBQPCSFQJVCyBxUCRfGQWFA2sv1QlOFEhUBRPGiMQMlGVAgUxYMjAEEVAiWLUCFQFA+aPqpkIyoEiuIhoAvKlCxEhUBRPKgAKNmICoGieHBLTOidoWQRerkrigdB1xEo2YcKgaJ4cCwBFQIlm1AhUBQPmjWkZCMqBIrig+qAkk2oECiKB2NvgqquISWbyEgIRKRURBaLyEb7d0madk+LSIOI/C3l+BQReUNENonIIyISyaQ/ipIpcVsJVAiUbCJTi+BGYKkxZjqwlPSbz/8YuNrn+A+BnxtjjgT2A9dk2B9FyQhHCFQHlGwiUyGYBzxgP34AuNSvkTFmKdDsPSZWUZePDGOHKwAABlFJREFUAI/1935FGS7UNaRkI5kKQbkxptp+vAcoP4D3jgYajDEx+/lOYGKG/VGUjHCFQKNnShbR7+b1IrIEGOfz0s3eJ8YYIyJmsDrm04/rgOsAJk+ePFRfo2Q5rmtI84aULKJfITDGnJ/uNRHZKyLjjTHVIjIe2HcA310HFItIyLYKKoBdffRjAbAAoKqqasgER8luHJdQbiQ4wj1RlOEjUwN4ITDffjwf+MtA32iMMcBzwJUf5P2KMhQcO7GIf/3Ikdx11Ukj3RVFGTYyFYI7gQtEZCNwvv0cEakSkXudRiLyEvBH4DwR2SkiF9ovfRu4QUQ2YcUM7suwP4qSESLCDXOPYtyonJHuiqIMG/26hvrCGFMHnOdzfDlwref5WWnevwU4JZM+KIqiKJmhuRGKoihZjgqBoihKlqNCoCiKkuWoECiKomQ5KgSKoihZjgqBoihKlqNCoCiKkuWIMYdetQYRqQG2f8C3jwFqB7E7g4X268DQfh0Y2q8D43Dt1xHGmLLUg4ekEGSCiCw3xlSNdD9S0X4dGNqvA0P7dWBkW7/UNaQoipLlqBAoiqJkOdkoBAtGugNp0H4dGNqvA0P7dWBkVb+yLkagKIqiJJONFoGiKIri4bAVAhH5sYi8JyJrROTPIlLsee0mEdkkIus9eyMgIhfZxzaJyI1D1K9PiMg6EYmLSJXneKWItIvIW/bPPZ7X5ojIWrtfd4kM/s7q6fplvzZi58unn7eKyC7Pefpof/0cLkbifPTRl232NfOWiCy3j5WKyGIR2Wj/LhmGftwvIvtE5G3PMd9+iMVd9vlbIyKzh7lfI35ticgkEXlORN6x78ev2seH9pwZYw7LH2AuELIf/xD4of14JrAaiAJTgM1A0P7ZDEwFInabmUPQr2OAo4DngSrP8Urg7TTveRM4DRDg78DFw9ivET1fPv28FfiGz3Hffg7j9TYi56OP/mwDxqQc+xFwo/34RueeGOJ+nA3M9l7b6foBfNS+vsW+3t8Y5n6N+LUFjAdm248LgQ329w/pOTtsLQJjzDPG2gsZ4HWsPZEB5gEPG2M6jTFbgU1Ym+OcAmwyxmwxxnQBD9ttB7tf7xpj1g+0vb0XdJEx5nVj/ed/D1w6jP0a0fN1AKTr53BxsJ0PP+YBD9iPH2AIrqNUjDEvAvUD7Mc84PfG4nWsPc3HD2O/0jFs15YxptoYs9J+3Ay8C0xkiM/ZYSsEKXwBSzXBOqk7PK/ttI+lOz6cTBGRVSLygog4u7pNtPsyUv06GM/X9bYZfL/HvTHS/7+R/v5UDPCMiKwQkevsY+XGmGr78R6gfGS6lrYfB8M5PGiuLRGpBE4C3mCIz1lGW1WONCKyBBjn89LNxpi/2G1uBmLA/x5M/fKhGphsjKkTkTnAEyIy6yDo17DTVz+B/wJuwxrobgN+iiX0SjJnGmN2ichYYLGIvOd90RhjRGTEUwYPln7YHDTXlogUAI8DXzPGNHnDgkNxzg5pITDGnN/X6yLyeeAfgPNstwrALmCSp1mFfYw+jg9qv9K8pxPotB+vEJHNwAy7DxWepsPaL4bhfKUy0H6KyH8Df7Of9tXP4WCkvz8JY8wu+/c+Efkzlitjr4iMN8ZU2+6DfSPUvXT9GNFzaIzZ6zweyWtLRMJYIvC/xpg/2YeH9Jwdtq4hEbkI+BZwiTGmzfPSQuAqEYmKyBRgOlYwdhkwXUSmiEgEuMpuO1z9LRORoP14qt2vLbY52CQip4k1LfgcMJyz94PqfKX4Py8DnKyPdP0cLkb0+vEiIvkiUug8xkqceNvuz3y72XyG9zrykq4fC4HP2ZkwpwGNHnfIkHMwXFv2PX4f8K4x5meel4b2nA1F5Ptg+MEK6OwA3rJ/7vG8djNW5H89ngwcrAj8Bvu1m4eoX5dh+fE6gb3AIvv4FcA6u68rgY973lOFdVFuBn6NvRBwOPo10ufLp58PAmuBNfZNML6/fg7jNTfs5yNNP6ZiZbmstq+pm+3jo4GlwEZgCVA6DH35Pyy3Z7d9fV2Trh9YmS+/sc/fWjzZa8PUrxG/toAzsVxTazxj10eH+pzpymJFUZQs57B1DSmKoigDQ4VAURQly1EhUBRFyXJUCBRFUbIcFQJFUZQsR4VAURQly1EhUBRFyXJUCBRFUbKc/w85ucS1yltGuQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with pyasdf.ASDFDataSet(cc_h5,mode='r') as ds:\n",
    "    data_type = ds.auxiliary_data.list()\n",
    "    path = ds.auxiliary_data[data_type[0]].list()\n",
    "    print(data_type,path)\n",
    "    \n",
    "    data = ds.auxiliary_data[data_type[0]][path[0]].data[:]\n",
    "    para = ds.auxiliary_data[data_type[0]][path[0]].parameters\n",
    "    \n",
    "    # plot the waveform again\n",
    "    plt.plot(tvec,data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The end.\n",
    "\n",
    "We hope you enjoy it! Most of the core steps of NoisePy are included here for illustration, and NoisePy is simply adding more loops for source and receiver stations and use some tricks (including embedding parallel functionality) to speed up the general performance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
